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I.  INTRODUCTION 

The  economics  associated  with  ship  operations  have 
necessitated  an  examination  of  the  losses  associated  with 
the  motion  of  an  automatically  steered  ship  in  a  seaway. 

Four  major  areas  where  fuel  losses  occur  during  the 
operation  of  a  ship  have  been  identified  [Ref.  1,  2].  These 
areas  on  existing  steam/diesel  tankers  are  shown  below: 

•  Power  plant  and  auxiliaries 

•  Propeller  efficiency 

•  Hull  resistance 

•  Steering  and  navigation 

An  optimized  autopilot  design  would  provide  effective 
steering  control  with  associated  cost  savings  due  to 
reducing  fuel  consumption. 

An  appropriate  computer  model  which  represents  the  ship 
is  necessary  for  studies  leading  to  appropriate  controller 
design.  Chapter  2  introduces  the  development  of  two  of  these 
models . 

Chapter  3  addresses  the  formulation  of  a  performance 
criterion  which  represents  the  added  drag  due  to  steering  of 
the  ship. 

Using  the  equations  of  motion  as  a  model  of  the  ship  and 
a  function  minimization  subroutine  we  proceed  to  the 
controller  design  for  regular  seas  (deterministic  model  for 
the  seaway)  in  Chapter  4,  and  for  irregular  seas  (nondeter- 
ministic  model)  in  Chapter  5.  The  function  minimization 
subroutine  used  was  BOXPLX  and  was  programmed  by  R.  R. 
Hilleary  of  the  Naval  Postgraduate  School  Computer  Center 
[Ref.  3].  It  will  find  the  minimum  of  any  arbitrary  func- 
tion, linear  or  nonlinear,  subject  to  explicit  constraints 
of  the  variables  or  implicit  constraints  on  functions  of  the 
variables . 
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Chapter  6  introduces  another  function  minimization 
subroutine  appropriate  for  onboard  use. 

An  adaptive  control,  which  updates  the  controller  param- 
eters when  the  environmental  conditions  or  the  ship's  course 
change,  is  studied  in  Chapter  7. 

Conclusions  drawn  from  these  experiments  and  recommenda- 
tions for  future  studies  are  addressed  in  Chapter  8. 
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II.  COMPUTER  MODELS  OF  THE  SHIP 

A  nontrivial  part  of  any  control  problem  is  modelling 
the  process.  Thus,  an  appropriate  computer  model  which 
represents  the  ship  is  necessary.  The  best  representation  of 
the  ship's  steering  dynamics  is  a  Taylor's  series  expansion 
of  the  force  and  moment  relationships  around  a  selected 
steady  state  operating  point.  The  equations  obtained  in  this 
way  are  known  as  the  equations  of  motion  [Ref .  4] ,  and  the 
formulation  in  the  computer  program  is  indicated  in  Appendix 
A.  This  computer  program  was  developed  by  using  known  avail- 
able data  for  the  SL-7  containership  and  by  implementation 
of  the  scheme  in  Figure  2.1  [Ref.  5]. 

In  this  scheme  the  function  minimization  subroutine  is 
fed  by  the  yaw  error  \y  and  rudder  angle  Q     ,       computes  the 

re 

performance   criterion  J   and   adjusts   the  controller   free 
parameters  in  order  to  minimize  J. 

A  second  model  for  the  ship- steering  dynamics  represen- 
tation is  the  Nomoto  model.  Figure  2.2  indicates  the  second 
and  third  order  Nomoto  transfer  functions  while  Figure  2.3 
indicates  the  appropriate  scheme  used  for  obtaining  these 
models  from  the  equations  of  motion.  Appendix  A  includes 
the  computer  program  used  for  the  Nomoto  third  order  model 
determination . 

A  yaw  command  is  applied  as  input  in  the  scheme  in 
Figure  2.3  and  the  difference  of  the  signals  \u  and  U^_~  is 
fed  to  the  function  minimization  subroutine  wnich  attempts 
to  adjust  the  free  parameters  of  the  Nomoto  plant  in  order 
to  minimize  the  performance  criterion  J. 

Simulation  runs  indicate  that  the  resulting  Nomoto 
models  are  obtained  with  resulting  J  close  to  zero. 
However  ,  in     this    study  .  the    equations    of    motion 
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representation  was  adopted  because  the  system  is  dynamic 
and  use  of  the  Nomoto  model  representation  implies  addi- 
tional computer  use.  On  the  other  hand,  frequency  domain 
studies  were  carried  out  using  the  Nomoto  representation 
since  this  representation  is  easier  handled. 
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III.  AN  ADEQUATE  PERFORMANCE  CRITERION 

A.   CRITERION  BASED  ON  TRUE  ADDED  RESISTANCE 

The  performance  criterion  which  characterizes  propulsion 
losses  due  to  steering  may  be  shown  to  be  that  derived  from 
excess  power  consumption  per  unit  distance  caused  due  to 
steering  [Ref.  1,  6].  The  added  resistance  due  to  steering 
can  be  related  to  the  surge  or  thrust  equation  where  the 
total  instantaneous  surge  relevant  to  steering  is 

AX=[m+(p/2)LAX^r]yr  +  l/2[(p/2)AX/vJv2  (3.1) 

.+1/2[(P/2)AX^  U2]52 

where   m  =  mass  of  ship 

P  -  density  of  sea  water 

L  =  ship's  length  between  perpendiculars 

A  =  L2 

U  =  ship's  water  speed 

v  =  sway  velocity 

r  =  yaw  rate  of  ship 

0  =  rudder  angle 

X'  =  force  coefficient  due  to  yaw/sway  (posi- 


tive ) 
Xij-  =  force  coefficient  due  to  rudder  angle 

(negative ) 
X'   =  force  coefficient  due  to  sway 


Since  the  sway  velocity  of  the  ship  is  small  we  can 
neglect  the  term  which  includes  the  square  of  the  sway 
velocity  in  the  previous  equation.  From  this  the  mean  surge 
relevant  to  steering  may  be  written  as 
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AX=[m+(  p/2)LAX^](uara/2)cos((^  -  if   )  (3.2) 

+  [(P/2)AX^U2](  52/2)      V         r 


Where      u   =  amplitude  of  sway  velocity 
r   =  amplitude  of  yaw  rate 


Qq    =  amplitude  of  rudder  angle 

(y  -  (fl  =    phase  difference  between  sway  and 
V      r 

yaw  rate 


A   performance  criterion   for   added   resistance  due   to 
steering  may  be  formulated  as 

J=lim  (l/2T)/(- J  vr+yu2  5  2)^t  (3.3) 

T—oo    Jq 

where  Q    and  'V  are  constants 


Accurate  knowledge  of  the  nonlinear  coefficients  X'  and 
Xsjt  is  required  for  the  accuracy  of  such  a  criterion.  In 
addition  the  criterion  itself  suffers  from  the  disadvantage 
that  sway  velocity  measurements  are  not  available. 
Normalizing  the  last  equation  the  performance  criterion  will 
be 

■T 

J    =lim  (l/2T)/(-Xvr+52)dt  (3.4) 

normT-*oo         Jq 

where  A={2[m+  (  p  /2)LA]X^  }  /  [  (ft  /2  )X^  U2] 

Table  I   indicates  the   values  of   A  for  the   operating 
range  of  speed  of  the  ship  studied. 


17 


TABLE  I 
Weighting  factor 


X 
X' 


Ship's  speed  (knots) 

16  21.5350 

23  10.4215 

32  5.3900 


B.   CRITERION  BASED  ON  APPROXIMATE  ADDED  RESISTANCE 

Empirical  criteria  based  on  an  approximation  to  added 
resistance  may  also  be  derived.  A  semiempirical  criterion 
for  measuring  the  relative  performance  was  developed 
[Ref.  7],  based  on  the  assumption  of  small  amplitude  oscil- 
lations around  the  steady-state  pivot  point  of  the  ship 
during  yawing  at  the  ship/ steering  system  natural  frequency. 
This  may  be  extended  and  an  alternative  criterion  for  added 
resistance  will  be 

T 

J=lim  (l/2T)|(y\i/^2+  52)dt  (3-5) 

7— co 

where     t 

A=A(J  =  [2m(l  +  X;r)(0P/L)CJ2]/[(P/2)LX^U2] 

V=  [(P/2)LAX^.]/m 

OP  =  distance  from  center  of  gravity  to  pivot 

center 
OJ    =    natural  frequency  (closed  loop  ship 

steering  control) 
yi   =  ship's  perturbation  yaw  angle 

The  values  of  A  as  a  function  of  ship's  speed  are  given 
by  Table  II. 

A  closed  loop  system  natural  frequency  (jj  of  around  0.05 
rads  per  sec   has  the  potential  to  attenuate   the  effects  of 
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TABLE  II 

/ 

Weighting  factor 

X 

nots  ) 

X 

6,720 
3,251 
1,680 

Ship's  speed  (knots) 

23 
32 


seaway  disturbance  in  the  range  of  encounter  angles  where 
added  resistance  due  to  steering  is  important  [Ref.  6].  The 
weighting  factor  for  the  operating  range  of  the  ship  is 
shown  in  Table  III. 


TABLE  III 
Weighting  factor 


X 


Ship's  speed  (knots) 


X 


6  16.796 

23  8.128 

32  4.2 


Equation  3.5  is  used  as  a  performance  criterion  for  this 
study.  It  is  an  approximation  but  it  is  convenient  for 
onboard  use  since  ship's,  perturbation  yaw  angle  \p  and 
rudder  angle  Q      are  measurable. 

C.   WEIGHTING  FACTOR  STUDY. 

The  weighting  factor  A  given  by  Table  III  used  in  equa- 
tion 3.5,  plays  an  important  role  in  terms  of  the  optimal 
controller  parameters  determination.  Some  investigation  is 
necessary  in  order  to  verify  the  accuracy  of  the  results, 
since  the  values  of  A  of  Table  III  are  determined  based  on 
the  assumption  that  the  closed  loop  system's  natural 
frequency  is  around  0.05  rads  per  sec  [Ref.  1].  Frequency 
domain  techniques  were  used  for  this  purpose.  Using  the 
Nomoto  third  order  model  representation  of  the  ship  and 
available  controller  parameters  from  Chapter  4  for  sea  state 
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4,  encounter  frequency  1.5  rads  per  sec,  encounter  angle 
150°  and  ship's  speed  23  knots  we  found  that  the  closed  loop 
bandwidth  of  the  system  is  0.04  rads  per  sec,  as  indicated 
in  Figure  3.1,  which  is  not  close  enough  to  0.05  rads  per 
sec . 

For  the  same  sea  conditions  and  ship  speed,  with  the 
assumption  that  the  closed  loop  natural  frequency  of  the 
system  is  not  0.05  rads  per  sec  but  0.04  rads  per  sec,  a  new 
value  A=5.734  was  obtained  and  the  frequency  domain  tech- 
niques result  in  a  new  bandwidth  0.035  rads  per  sec  as  is 
indicated  in  Figure  3.2. 

Clearly,  the  values  of  A  given  by  Table  III  and  used  in 
this  study  are  not  the  best.  Unfortunately,  since  the  full 
hydrodynamic  coefficients  of  the  SL-7  containership  are  not 
known  we  can't  develop  the  surge  equation  and  thus  it  is 
still  impossible  to  determine  accurate  values  for  the 
weighting  factor  A. 
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Figure  3.1    Closed  Loop  Bandwidth  for  A=8  .  12. 
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Figure  3.2    Closed  Loop  Bandwidth  for  A  =  5.734 
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IV.  REGULAR  SEAS  L    CONTROLLER  DESIGN 

We  have  already  defined  a  suitable  and  sufficiently 
accurate  ship  computer  model  and  the  system's  performance 
criterion,  as  well.  The  remaining  task  is  to  determine  a 
representation  of  the  external  disturbances  imparted  to  the 
ship  by  the  sea,  before  the  system's  performance  in  a  seaway 
can  be  evaluated.  A  correct  model  of  the  seaway  itself  is 
essential  to  representative  modeling  of  forces  and  moments 
exerted  on  the  ship  by  it. 

At  this  point  we  will  use  the  regular  sea  model  as  sea 
representation.  The  properties  of  regular  seas  are  well 
defined.  The  wave  crests  are  assumed  to  be  straight,  infi- 
nitely long,  parallel  and  equally  spaced  with  constant  wave 
height.  The  waves  progress  in  a  direction  perpendicular  to 
the  crest  line  at  a  uniform  velocity.  However  the  sea  is 
never  regular.  It  is  a  random  phenomenon  where  waves  are 
continually  changing  in  height,  length  and  breadth  [Ref.  8]. 

The  forces  exerted  by  the  regular  sea  have  the  form 

F  =  (j0R/cos(q?t+#/  )  (4.1) 

where   ^n~    significant  wave  height 
R:  =  exciting  force 
CJU=  encounter  frequency 
Xfj-    phase  angle 
The  correspondence  between  sea  state   and  wave  height  is 
indicated  in  Table  IV  [Ref.  9]. 

The  exciting  forces  R-  for  different  encounter  frequen- 
cies and  encounter  angles  were  obtained  from  the  sea  state 
generator  program  [Ref.  10], 
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TABLE  IV 
Sea  state  vs  Wave  height 


Sea  state 
(Beaufort  scale) 
0 
1 
2 
3 
4 
5 
6 
7 
8 
9 


Range  for  wave  height 
(Feet) 

0.0 

0.32 
0.65-0.98 
1.96-3.28 
3.28-4.92 
6.56-8.20 
9.84-13.1 
13.1-18.2 
18.2-24.6 
23.0-32.9 


Appendix  B  indicates  the  regular  seastate  formulation  in 

the  FORTRAN  program  used  for  obtaining  the  controller  param- 
eters . 

The  controller  used  in  the  entire  study  has  one  pole-one 

zero   and   the   form   is  indicated   in   Figure   4.1.    This 

controller  seems  to  have  the  best  performance  in  calm  waters 
and  in  seaway  [Ref.  5], 


% 

„ 

8 

K/0+T/s) 
(!+T2s) 

Figure  4.1    Controller  Used  in  this  Study 

The  optimized  controller  parameters  and  the  cost  J  for 
23  knots  speed,  sea  states  4-6-7-9,  different  encounter 
angles  and  various  encounter  frequencies  are  indicated  in 
Tables  V, VI, VII  and  VIII. 
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Studying  the  Tables  V  through  VIII  we  can  draw  the 
following  conclusions: 

•  For  a  particular  encounter  angle  and  encounter 
frequency  the  higher  the  sea  state  the  higher  the  cost. 

•  For  the  same  sea  state  the  cost  becomes  smaller  for 
higher  encounter  frequencies . 

•  For  encounter  frequency  0.2  rads  per  sec  the  maximum 
cost  occurs  at  60°  encounter  angle  for  all  tested  sea 
states . 

•  For  0.6  and  0.75  rads  per  sec  encounter  frequency  the 
maximum  cost  occurs  at  120°  encounter  angle  for  all 
tested  sea  states . 

•  For  1.5  rads  per  sec  encounter  frequency  the  maximum 
cost  occurs  at  90°  encounter  angle  for  all  tested  sea 
states . 

•  For  0.4  rads  per  sec  encounter  frequency  the  maximum 
cost  occurs  at  90°  encounter  angle  for  sea  states  4,6 
and  7  while  at  sea  state  9  the  maximum  cost  occurs  at 
60°  encounter  angle. 

Appendix  C  provides  the  computer  program  necessary  to 
achieve  the  system's  response.  Some  typical  responses  are 
indicated  in  Figures  4.2,  4.3,  4.4,  4.5,  4.6,  4.7,  4.8,  4.9. 
It  is  obvious  that  as  the  sea  state  goes  to  heavier  seas  the 
rudder  and  yaw  perturbations  become  larger. 

An  attempt  to  determine  how  accurate  the  controller 
parameters  must  be  for  a  particular  situation,  leads  to  the 
conclusion  that  high  accuracy  isn't  required.  Keeping  two 
parameters  fixed  each  time  and  vary  the  third  we  can  see 
(Figures  4.10,  4.11,  4.12,  4.13,  4.14)  that  the  cost  doesn't 
change  appreciably  in  the  vicinity  of  the  actual  value. 

Figures  4.2,  4.3,  4.4,  4.5,  4.6,  4.7,  4.8,  4.9  indicate 
that  the  yaw  and  rudder  excursions  are  less  than  1°.  This 
just  seems  strange,  though  it  may  be  because  of  the  opti- 
mization of   the  filter.    We  tried   to  investigate   that  by 
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using  the  optimal  filter  for  sea  state  9,  encounter 
frequency  1.5  rads  per  sec,  encounter  angle  120°  and  run  it 
in  sea  state  4,  keeping  the  same  encounter  frequency  and 
angle.  The  yaw  and  rudder  excursions,  even  if  they  became 
larger,  remained  less  than  1°.  The  parameters  of  those  two 
filters  are  close  and  the  reason  might  be  the  flatness  of 
the  cost  surface.  Second  attempt  led  to  more  interesting 
results.  Using  the  same  filter  and  run  it  in  sea  state  4, 
encounter  frequency  0.4  rads  per  sec  and  encounter  angle 
060°,  the  system  becomes  unstable  (Figures  4.15,  4.16). 
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TABLE  V 

Optimal  Controller  Parameters  for  Regular  Sea 

Sea  State  4 


Encounter  Frequency  0.2  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees ) 

0        0.5221067  66.3312231  12.8332741  0.609E-33 

30         1.0488701  61.9309387  15.9266357  2.4819 

60         1.2036362  54.5533295  16.0245972  3.612223 

90         1.3178062  49.7426453  14.8329315  2.703524 

120         1.3984699  46.9797058  13.9525757  1.355578 

150         1.4502153  45.4263306  13.3351599  0.3567196 

180         0.7195223  25.2119598  14.1219782  0.345E-28 

Encounter  Frequency  0.4  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees ) 

0         0.5221067  66.3312231  12.8332741  0.517E-35 

30         0.7234516  74.7846533  47.9879893  0.54673 

60         0.8730211  92.8420868  50.0454407  1.194746 

90         2.8910360  28.9679871  10.0176315  2.536161 

120         2.6232796  27.9702454  8.8327761  0.2349457 

150         2.6408129  27.7937927  8.6838455  0.0545772 

180         0.7195223  25.2119598  14.1219782  0.536E-32 

Encounter  Frequency  0.6  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.453E-31 

30         2.0506487     1.2107763  5.5998001  0.0028366 

60         1.0504951     0.2329917  19.0010876  0.0031525 

90         2.1496201     1.2607498  5.3318434  0.0790777 

120         2.1221027     0.9715905  6.9064713  0.0796856 

150         1.9617786     0.8480816  8.3953667  0.0341976 

180         0.7195223  25.2119598  14.1219782  0.147E-39 

Encounter  Frequency  0.75  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.711E-35 

30         2.4342451     0.8039263  7.4397717  0.0017128 

60         2.3829517     0.8094406  8.6272535  0.0042339 

90         2.0794163     0.4690621  16.9594727  0.0442135 

120         1.9938784     0.2545378  22.4892426  0.1164415 

150         1.9387684     0.2628353  23.8461609  0.0379958 

180         0.7195223  25.2119598  14.1219782  0.312E-23 

Encounter  Frequency  1.5  rads  per  sec 

encounter        Kl           Tl  T2             J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.143E-35 

30         1.8784037     0.6553955  5.4344263  0.0000847 

60         2.4894753     0.5992758  5.1160698  0.0014724 

90         2.5899792     0.6663168  3.1813316  0.0241299 

120         1.8784037     0.5395500  10.4344263  0.0028635 

150         2.4478331     0.5559916  5.9208412  0.0015984 

180         0.7195223  25.2119598  14.1219782  0.451E-28 
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TABLE  VI 

Optimal  Controller  Parameters  for  Regular  Sea 

Sea  State  6 


Encounter  Frequency  0.2  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees ) 

0        0.5221067  66.3312231  12.8332741  0.322E-33 

30         1.0966187  62.3766327  19.3941193  8.622314 

60         1.2665367  55.1072540  19.2263336  11.84099 

90         1.3603411  49.3086395  16.7315979  9.270265 

120        1.4170542  46.4313202  14.7203903  5.004774 

150         1.4533644  45.2830963  13.6198730  1.396385 

180        0.7195223  25.2119598  14.1219782  0.678E-28 

Encounter  Frequency  0.4  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees ) 

0         0.5221067  66.3312231  12.8332741  0.122E-31 

30         0.7536127  12.4472111  4.8151121  2.3798 

60         1.1621914     1.7206783  5.6173820  4.143427 

90         2.8681517  28.5688019  11.3736725  7.7779746 

120         2.6123600  28.0108032  9.1387405  0.8634676 

150         2.6361399  27.8540497  8.7665482  0.2140626 

180         0.7195223  25.2119598  14.1219782  0.345E-25 

Encounter  Frequency  0.6  rads  per  sec 

encounter        Kl           Tl  T2           J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.806E-35 

30         2.0449467     1.2405663  5.5909785  0.0113364 

60         1.0594482     0.1432046  18.5202637  0.0126108 

90         2.1537333     1.1794500  5.9987049  0.2766824 

120         2.1456242     0.9514294  7.3327799  0.3125796 

150         1.9752455     0.8255181  8.5351868  0.1364259 

180         0.7195223  25.2119598  14.1219782  0.112E-23 

Encounter  Frequency  0.75  rads  per  sec 

encounter       Kl           Tl  T2           J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.753E-32 

30         2.4413719     0.7599964  7.3472633  0.0068484 

60         2.4142313     0.7818718  8.5500679  0.0169234 

90         2.0592680     0.3840635  19.2782745  0.1523162 

120         2.0695038     0.2843146  22.5692444  0.4640285 

150         1.9513931     1.2351496  23.5140228  0.1518951 

180         0.7195223  25.2119598  14.1219782  0.691E-28 

Encounter  Frequency  1.5  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.321E-33 

30         1.7606564     0.4668925  12.3475494  0.0003390 

60         2.5002985     0.5747030  5.3262844  0.0058868 

90         2.6015043     0.6543975  3.2408133  0.0946725 

120         2.1809998     0.4698086  10.9160814  0.0114672 

150         2.4302263     0.5666089  6.0247307  0.0063947 

180         0.7195223  25.2119598  14.1219782  0.344E-23 
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TABLE  VII 

Optimal  Controller  Parameters  for  Regular  Sea 

Sea  State  7 


Encounter  Frequency  0.2  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees ) 

0        0.5221067  66.3312231  12.8332741  0.911E-35 

30         1.1839037  68.3722687  27.3912964  19.62723 

60         1.3688898  65.5752258  30.9693604  24.53816 

90         1.4492817  53.9667511  23.6652069  20.31325 

120         1.4652090  46.2668457  17.1322327-  12.58828 

150        1.4653692  44.8378906  14.1106033  4.036665 

180         0.7195223  25.2119598  14.1219782  0.912E-32 

Encounter  Frequency  0.4  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.341E-30 

30         0.6452138  38.7243542  11.6571761  0.0524352 

60         2.2381306     1.0741718  9.7088852  10.22498 

90         2.7780743  30.6558838  16.9282227  11.95129 

120         2.5894566  28.1353302  10.0156784  2.098495 

150         2.6306906  27.8760681  8.9550962  0.6211755 

180         0.7195223  25.2119598  14.1219782  0.176E-35 

Encounter  Frequency  0.6  rads  per  sec 

encounter        Kl          Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.157E-31 

30         2.0501146     1.2125063  5.5930214  0.0346353 

60         1.0624285     0.1251755  18.0873718  0.0386276 

90         2.1488247     0.9770966  8.1855469  0.6328694 

120         2.1970577     0.8763103  8.4206886  0.9178923 

150         2.0124264     0.8064904  8.9983368  0.4150548 

180         0.7195223  25.2119598  14.1219782  0.441E-28 

Encounter  Frequency  0.75  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees ) 

0         0.5221067  66.3312231  12.8332741  0.238E-28 

30         2.4456072     0.8701959  7.6628313  0.0209509 

60         2.4188919     0.8008499  8.6116581  0.0517269 

90         1.9072247     0.1301596  28.8665771  0.4958954 

120         2.2209492     0.3453745  22.8497772  1.398130 

150         2.0318069     0.2577103  23.6702576  0.4641950 

180         0.7195223  25.2119598  14.1219782  0.542E-21 

Encounter  Frequency  1.5  rads  per  sec 

encounter        Kl           Tl  T2             J 
angle (degrees ) 

0         0.5221067  66.3312231  12.8332741  0._712E-32 

30         2.0834208     0.4371362  15.1762695  0.0010384 

60         2.4635468     0.5746776  5.3325920  0.0180051 

90         2.6105270     0.6474890  3.4996614  0.2765892 

120         2.1780138     0.4588630  11.4343033  0.0352355 

150         2.4314880     0.5580992  6.0874767  0.0195939 

180         0.7195223  25.2119598  14.1219782  0.413E-26 
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TABLE  VIII 

Optimal  Controller  Parameters  for  Regular  Sea 

Sea  State  9 


Encounter  Frequency  0.2  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees ) 

0        0.5221067  66.3312231  12.8332741  0.493E-31 

30         1.2344990  97.2906342  53.3980713  31.88802 

60         1.5019569  77.3820190  45.6973572  30.87321 

90         1.5020103  77.2400513  45.6112671  30.87321 

120         1.5420017  51.2350922  23.8943634  21.75188 

150         1.4879055  44.4281464  15.2255249  8.599576 

180        0.7195223  25.2119598  14.1219782  0.810E-39 

Encounter  Frequency  0.4  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees ) 

0         0.5221067  66.3312231  12.8332741  0.743E-33 

30         0.1618259  95.8238471  18.3397064  8.103665 

60         2.3103380  0.9894915  12.5459442  20.30595 

90         3.3482129  2.6129665  4.7000313  7.494904 

120         2.5573978  28.7854462  12.1667938  3.113852 

150         2.6185656  27.9137726  9.3553696  1.324692 

180         0.7195223  25.2119598  14.1219782  0.534E-33 

Encounter  Frequency  0.6  rads  per  sec 

encounter       Kl           Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.478E-31 

30         2.0493793  1.2093735  5.6404839  0.0820505 

60         1.1003008  0.1920759  18.7369995  0.0919840 

90         2.0454016  0.5407953  15.4096680  1.031631 

120         2.2776527  •   0.7766371  10.4350891  2.069717 

150         2.0829115  0.7725539  9.8721237  0.9771649 

180         0.7195223  25.2119598  14.1219782  0.117E-29 

Encounter  Frequency  0.75  rads  per  sec 

encounter        Kl           Tl  T2            J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.872E-35 

30         2.4433956  0.8601446  7.7009745  0.0497640 

60         2.4255047  0.7809458  8.7044067  0.1226490 

90         2.2933350  0.4751374  16.7149658  0.7522082 

120         2.3928595  0.4103206  22.9515076  3.146294 

150         2.1456175  0.2891768  23.5924225  1.097565 

180         0.7195223  25.2119598  14.1219782  0.611E-26 

Encounter  Frequency  1.5  rads  per  sec 

encounter       Kl           Tl  T2             J 
angle (degrees) 

0         0.5221067  66.3312231  12.8332741  0.242E-33 

30         2.0808840  0.4744592  16.3423157  0.0024732 

60         2.4694090  0.5775681  5.4030972  0.0042757 

90         2.6412249  0.6233797  3.9802380  0.6127556 

120         2.1955976  0.4414446  12.1491365  0.0084519 

150         2.4388838  0.5490999  6.1774960  0.0046688 

180         0.7195223  25.2119598  14.1219782  0.420E-30 
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V.  IRREGULAR  SEAS  2  CONTROLLER  DESIGN 

The  major  characteristic  of  the  sea  is  its  irregularity. 
This  irregularity  can  be  described  by  statistical  methods  by 
assuming  that  a  large  number  of  regular  (sinusoidal)  waves 
having  different  wavelengths,,  directions,  phases  and  ampli- 
tudes are  superimposed  to  form  the  randomly  varying  sea. 

The  presence  of  the  irregular  sea  was  obtained  by 
coupling  a  sea  state  generator  program  to  the  FORTRAN 
program  as  is  indicated  in  Appendix  D.  The  sea  state  gener- 
ator program  generates  added  mass  and  added  inertia  values 
as  function  of  the  encounter  frequency  and  also  calculates 
forces  and  moments  imparted  to  the  shiphull  by  the  sea.  The 
forces  and  moments  are  stored  in  a  look  up  table  which  was 
coupled  to  the  equations  of  motion.  The  irregular  sea  waves 
impinging  on  the  ship  contain  the  total  energy  density  spec- 
trum composed  of  many  frequencies  and  the  ship  responds  to 
an  average  value  of  added  mass  and  added  inertia,  while  in 
the  regular  sea  the  added  mass  and  added  inertia  were  known 
for  a  given  encounter  frequency.  We  decided  to  use  values 
for  added  mass  and  added  inertia  corresponded  to  encounter 
frequency  0.75  rads  per  sec,  since  the  energy  density  is 
maximum  in  the  vicinity  of  this  frequency  [Ref.  5].  This 
frequency  gave  us  values  representative  of  an  average  value 
for  added  mass  and  added  inertia. 

The  controller  used  for  this  study  was  the  controller 
described  in  Chapter  4  (Figure  4.1).  The  optimized 
controller  parameters  and  the  cost  J  for  sea  states  4,  6,  7, 
9  and  0°,-  30°,  60°,  90°,  120°,  150°,  180°  encounter  angles 
are  indicated  in  Table  IX. 

Studying  the  Table  IX  we  can  draw  the  following 
conclusions : 
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•  For  sea  states  6,  7,  9,  the  higher  the  sea  state  the 
higher  the  cost,  for  every  particular  encounter  angle. 

•  Comparing  costs  for  sea  states  4  and  6  we  discover  some 
anomaly.  The  cost  for  a  specific  encounter  angle  in  sea 
state  4  is  higher  than  the  cost  for  the  same  encounter 
angle  in  sea  state  6.  Logically,  we  expect  higher  cost 
for  higher  sea  state. 

•  The  reason  for  this  anomaly  may  be  the  method  we  used 
in  order  to  obtain  the  added  mass  and  added  inertia 
values.  The  average,  we  consider,  might  not  represent 
the  actual  average. 

Appendix  E  provides  the  computer  program  necessary  to 
achieve  the  system's  response.  Some  typical  responses  are 
indicated  in  Figures  5.1,  5.2,  5.3,  5.4. 
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TABLE  IX 
Optimal  Controller  Parameters  for  Random  Sea 

Sea  State  4 

encounter        Kl  Tl  T2 

angle (degrees) 


0 

30 

60 

90 

120 

150 

180 

0. 
1. 
0. 
0. 

0. 
0. 
0. 

6021814 
5121580 
,6298036 
,6452737 
,7485995 
,9101038 
,6021814 

60 
89 
79 
82 
85 
92 
60, 

.30849 
.85324 
.07199 
,68692 
.77544 
,71379 
,30849 

10, 
19, 
10, 
10, 
12, 
15, 
10, 

,02579 
,85960 
.29221 
,79342 
.37746 
,21078 
.02579 

0, 
0, 
0, 
0, 
0, 
0, 
0, 

,342E-24 

,  10719 

.054196 

,076713 

,137624 

,012319 

,189E-28 

Sea 

State  6 

encounter 

Kl 

Tl 

T2 

J 

angle (degrees 

) 

0 

30 

60 

90 

120 

150 

180 

0, 

1, 

0, 
0, 
2, 
0, 
0, 

,6021814 
,8743490 
,8662014 
.7370305 
,6737600 
.4874309 
.6021814 

60 
61 
90 
84 
138 
77 
60 

.30849 
.82320 
.72922 
.08502 
.06650 
.86977 
.30849 

10, 
32 

14, 

12 

48 

41 

10 

.02579 
.22498 
.44058 
.17457 
.52447 
.73848 
.02579 

0 

0 

0, 

0 

0 

0 

0 

.172E-34 

.044758 

.028238 

.018541 

,065028 

.012478 

.767E-23 

Sea  State  7 

encounter        Kl  Tl  T2           J 
angle(degrees ) 

0         0.6021814  60.30849  10.02579  0.142E-34 

30         1.7232471  66.44597  27.32489  0.41237 

60         1.8508301  91.39410  36.91092  0.377894 

90         3.7642412  62.27244  86.58169  0.258193 

120         1.8482047  99.64923  91.28040  0.225806 

150         0.8519831  67.02328  64.96774  0.037724 

L80         0.6021814  60.30849  10.02579  0.811E-21 

Sea  State  9 

encounter       Kl  Tl  T2            J 
angle (degrees) 

0         0.6021814  60.30849  10.02579      0.711E-35 

30         3.1908160  138.56101  71.44171      1.306741 

60         3.0888780  152.01630  72.66160      0.961709 

90         3.2440700  121.13566  105.98480      0.523769 

120         1.5461040  111.49130  99.64659      0.2101076 

150         0.3758357  73.88629  35.17305      0.5007348 

180         0.6021814  60.30849  10.02579      0.121E-31 
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VI.  MINIMIZATION  SUBROUTINE  FOR  ONBOARD  USE 

A.   GENERAL 

As  is  mentioned  earlier  in  Chapter  1  the  function  mini- 
mization subroutine  used  for  these  studies  was  BOXPLX .  This 
subroutine  will  find  the  minimum  of  any  function,  linear  or 
nonlinear,  subject  to  explicit  constraints  of  the  variables 
or  implicit  constraints  on  functions  of  the  variables.  It 
will  handle  a  maximum  of  25  variables  but  can  handle  up  to 
50  variables  with  user  modification. 

The  variables  in  BOXPLX  are  allowed  to  move  within  a 
feasible  region  (n-dimensional  space,  where  n  is  the  number 
of  variables)  defined  by  upper  and  lower  bounds  on  their 
values .  The  choices  -  for  upper  and  lower  bounds  for  the 
parameters  are  based  on  an  understanding  of  the  function  of 
each  coefficient  of  the  system.  Experience  indicates  that 
while  accurate  selections  of  these  bounds  are  not  necessary, 
intelligent  selection  of  these  as  well  as  the  starting 
values  (guesses)  can  considerably  reduce  the  computer  number 
of  trials  needed  for  solution  convergence.  This  conclusion 
was  drawn  trying  to  obtain  the  controller  parameters  for 
Tables  V,  VI,  VII,  VIII  in  Chapter  4.  The  function  minimiza- 
tion subroutine,  when  starting  the  minimization  process  with 
arbitrary  chosen  guesses,  required  more  than  100  trials  for 
convergence  while  by  choosing  guesses  close  to  the  optimal 
parameters  required  more  than  50  and  less  than  100  trials. 
Considering  that  every  trial  lasted  600  seconds  (10  minutes) 
and  the  function  minimization  subroutine  requires  about  60 
samples  (trials)  before  telling  us  it  had  found  the  minimum 
this  would  mean  10  hours  for  the  control  to  adjust  itself. 
For  obvious  reasons,  such  operation  is  not  acceptable  for  on 
board  use. 
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For  obvious  reasons,  such  operation  is  not  acceptable  for  on 
board  use. 

B.  ATTACKING  THE  PROBLEM 

We  started  to  investigate  ways  to  improve  this.  These 
efforts  include: 

•  Finding   a    more   efficient    function   minimization 
subroutine 

•  Studying  the  flatness  of  the  cost  surface 

•  Reducing  sampling  time 

C.  SOLVING  THE  PROBLEM  - 

Switching  to  another  function  minimization  subroutine  we 
found  that  the  new  one  (ZXMWD)  suffered  from  the  same  disad- 
vantages . 

The  experiments  carried  out,  more  than  two  hundred, 
indicated  that  the  cost  surface  is  really  flat.  The  BOXPLX 
after  a  few  trials  started  to  focus  on  the  minimum  but 
before  it  converged,  it  needed  more  than  50  trials,  even  if 
the  guesses  were  close  to  the  optimal.  The  reason  is  the  way 
BOXPLX  itself  tries  to  find  the  minimum  of  a  function  of  NV 
variables.  It  converges  when  the  cost  FE  remains  unchanged 
for  2*NV  consecutive  trials  with  accuracy  10s.  An  effort  to 
modify  this  termination  criterion  in  terms  of  the  consecu- 
tive trials  was  successful.  Table  X  indicates  the  comparison 
between  modified  and  unmodified  BOXPLX,  for  sea  state  6, 
encounter  frequency  1.5  rads  per  sec  and  encounter  angle 
120°.  As  we  can  observe  in  Table  X  the  cost  in  each  case 
remains  almost  the  same  while  the  trials  required  for 
convergence  are  dependent  on  the  guesses  made  and  the  termi- 
nation criterion  established. 

The  value  of  the  cost  is  in  general  the  summation  of 
incremental  contributions   for  each  integration  step   and  is 
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TABLE  X 
First  Modification  in  BOXPLX 


BOXPLX 


Unmodified 

Unmodified 

Modified 

Modified 

Modified 

Modified 


Guesses 

Trials 

Arbitrary 

100 

Optimal 

72 

Arbitrary 

42 

Optimal 

11 

Arbitrary 

7 

Optimal 

1 

Terminat  ion 
Criterion 

6 
6 
3 
3 
2 
.2 


Cost 


0.01146720 
0.01146720 
0.01146812 
0.01146720 
0.01157926 
0.01146720 


therefore  dependent  on  the  total  time  of  the  simulation. 
This  is  important  in  that  the  optimal  gain  coefficients 
arrived  at  in  this  manner  are  not  optimal  for  steady-state 
performance  but  only  for  the  time  frame  covered.  This  should 
be  adequate  provided  the  time  frame  selected  is  long  rela- 
tive to  the  time  required  for  the  initial  condition  response 
to  die  out.  This  is  the  reason  Reid  has  chosen  time  frame 
600  seconds  [Ref.  1,2].  Of  course  this  time  period  is  large 
and  we  expect  steady-state  behaviour  faster  than  10  minutes. 
Simulation  studies  indicate  that  the  ship,  controlled  by  the 
controller  described  in  Figure  4.1,  reaches  the  steady-state 
situation  in  less  than  100  seconds.  So,  we  can  reduce  the 
600  seconds  time  frame  to  200  seconds,  safely.  This  is  very 
important  since  now  the  modified  function  minimization 
subroutine  BOXPLX,  uses  samples  of  200  seconds  long  instead 
of  600  seconds,  converges  in  less  than  30  minutes  which  is 
reasonable  for  on  board  use.  Since  the  value  of  the  cost  is 
dependent  on  the  time  frame  taken  we  expect  reduced  cost  for 
200  seconds  samples  long,  in  any  case. 

As  we  discussed  earlier  the  BOXPLX  compares  the  consecu- 
tive trials  with  accuracy  10s .  Since  we  do  not  need  so  big 
accuracy  a  second  modification  is  necessary.  We  decided  to 
change  the  existing  accuracy  to  lO**1*  . 

Table  XI  indicates  the  trials  required  for  BOXPLX 
convergence,  after  the  second  and  last  modification,  for  sea 
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TABLE  XI 
Second  Modification  in  BOXPLX 

Guesses      Trials      Termination  Cost 

criterion 

Arbitrary      17             3  0.003999837 

Optimal         8             3  0.003995277 

Arbitrary       5             2  0.004024354 

Optimal         1             2  0.003995247 


state  6,  encounter  frequency  1.5  rads  per  sec  and  encounter 
angle  120°.  By  comparing  Tables  X  and  XI  we  can  see  the 
further  improvement  for  the  function  minimization  subroutine 
convergence.  The  difference  in  the  cost  is  due  to  the 
different  time  frames  used.  The  modified  function  minimiza- 
tion subroutine  BOXPLX  is  indicated  in  Appendix  F. 
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VII.  ADAPTIVE  CONTROL 

A.  NECESSITY  OF  ADAPTIVITY 

The  plant,  the  system  which  is  supposed  to  be 
controlled,  is  normally  exposed- to  a  time  varying  environ- 
ment. The  ship's  speed,  the  wave  encounter  angle  and  the  sea 
state  are  changing  drastically  during  the  seaway.  Since 
changes  encountered  are  not  completely  predictable  an 
optimum  preprogrammed  time-varying  controller  is  not 
possible . 

If  we  assume-and  it  is  apparent  from  the  previous 
discussion  in  this  study-that  no  feasible  fixed  parameter 
controller  provides  acceptable  response  over  the  entire 
performance  spectrum,  it  is  necessary  that  some  means  is 
provided  for  adjusting  controller  parameters  according  to 
the  sea  conditions  and  ship's  operational  characteristics. 
Adaptive  control  is  thus  an  effort  to  extend  basic  optimum 
control  concepts  to  these  studies. 

B.  CANDIDATE  ADAPTIVE  SCHEMES 

Since  we  are  looking  for  1%  or  2%  savings  in  fuel  cost, 
we  have  to  feed  the  system  with  precise  information.  Exact 
knowledge  of  the  sea  state,  the  wave  encounter  angle  and  the 
ship's  ground  speed  is  vital  for  this  purpose.  Currently  the 
Navy  is  involved  in  a  program  that  will  provide  precision 
navigation  data.  Garcia  [Ref.  5],  provides  very  good  infor- 
mation on  this  subject. 

An  adaptive  control  scheme  is  indicated  in  Figure  7.1. 
Once  the  adaptive  part  of  the  scheme  is  set  up  the  sea 
state,  encounter  wave  angle  and  ship's  speed  are  fed  to  the 
filters  box   by  the   appropriate  sensors.    The  filters   box 
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includes  predetermined  optimal  filter  sets  for  different 
discrete  sea  states,  wave  angles  and  ship's  speed.  Actually, 
it  is  a  look  up  table  similar  to  those  of  Tables  V,  VI, VII 
and  VIII.  The  output  of  the  filter's  box  is  a  filter  which 
corresponds  to  discrete  conditions  close  to  those  fed  by  the 
sensors.  The  function  minimization  subroutine  accepting  the 
rudder  angle,  yaw  error  and  the  predetermined  filter  set  as 
initial  guesses  tries  to  obtain  the  optimal  filter  for  the 
exact  sea  state,  wave  encounter  angle  and  ship's  speed.  At 
the  same  time  the  plant  is  controlled  by  the  controller  with 
the  optimal  predetermined  set  of  parameters  for  conditions 
close  to  the  actual.  When  the  function  minimization  subrou- 
tine reaches  the  minimum,  it  supplies  the  controller  with  a 
new  set  of  parameters  which  is  the  optimal  set  for  the 
present  conditions. 

But,  what  happens  if  either  some  or  all  the  sensors 
provide  new  inputs  to  the  system?  A  decision  device  placed 
after  the  sensors  decides  whether  or  not  the  change  is 
appreciable.  This  device  compares  the  current  conditions 
with  those  used  to  obtain  the  controller  which  currently 
governs  the  system.  If  the  change  is  higher  than  some 
desired  percentage  then  a  new  predetermined  filter  is  passed 
to  the  controller  and  the  function  minimization  subroutine 
tries  to  find  the  'optimal  controller  parameters  for  the  new 
situation. 

From  the  scheme  of  Figure  7.1  we  can  eliminate  the 
predetermined  filters  device  which  provides  a  less  expensive 
system.  In  this  case  the  function  minimization  subroutine 
will  need  a  much  longer  time  to  determine  the  optimal  filter 
for  rapid  changes  in  course  and  speed,  even  if  we  assume 
that  rapid  changes  in  sea  state  don't  occur.  Finally,  the 
function  minimization  subroutine  might  work  continually  "on 
line"  as  is  indicated  in  Figure  7.2.  In  this  scheme  a  new 
controller  set  is  obtained  in  every  subroutine  iteration  and 
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some  constraints  for  the  controller  parameters  are  necessary 
in  order  to  avoid  operation  in  unstable  regions.  Thus,  the 
filter  supplied  by  the  subroutine  for  controlling  the  plant 
must  be  tested  for  characteristic  equation  roots  in  the 
right  half  S-plane.  In  that  case  a  default  option  must 
exist,  and  a  special  device  for  such  purposes  is  necessary. 
This  device  will  permit  change  in  the  filter  parameters  only 
when  the  new  set  still  keeps  the  system  in  a  stable  situ- 
ation. 

Whichever  adaptive  scheme  we  adopt,  we  must  provide  for 
manual  operation  for  the  system.  Manual  operation  will  be 
desired  in  the  following  cases: 

•  Arriving  ports 

•  Leaving  ports 

•  Restricted  waters 

•  Avoid  collision  in  open  seas 

•  Computer  down 

For  the  first  two  situations,  since  we  usually  expect  no 
heavy  seas,  the  optimal  filter  for  sea  state  1  seems  to  be 
more  suitable.  Also,  this  filter  can  serve  as  initial  condi- 
tion for  the  adaptive  schemes  described  before,  when  we  are 
leaving  ports.  For  the  rest  of  the  situations  the  last 
operating  optimal  filter  is  the  most  appropriate. 
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VIII.  CONCLUSIONS  AND  PvECOMMENDATIONS 

A.   CONCLUSIONS 

The  principal  conclusions  from  this  study  of  the  SL-7 
containership  as  they  related  to  steering  may  be  stated  as 
follows : 

•  It  is  evident  that  a  control  system  which  provides  the 
ship  heading  and  simultaneously  reduces  the  propulsive 
losses  does  exist  and  therefore  such  a  controller  saves 
fuel.  We  can't  conclude  how  much  the  savings  are, 
since  there  is  no  reference  for  comparison  between-  the 
conventional  autopilot  and  the  autopilot  which,  in 
addition,  holds  the  potential  for  reducing  the  propul- 
sive losses.  The  literature  says  that  savings  1%  or  2% 
is  possible. 

•  An  adaptive  controller,  that  minimizes  propulsion 
losses  as  ship  characteristics  and  .environmental  condi- 
tions change,  may  be  designed  using  a  self -  opt imalizing 
technique  employing  a  suitable  performance  criterion. 

•  Studying  every  particular  situation  we  conclude  that 
the  cost  surface  is  flat  and  therefore  accurate  deter- 
mination of  the  controller  parameters  is  not  required. 
So,  if  we  decide  to  use  the  adaptive  control  scheme  of 
Figure  7.1  we  don't  have  to  store  every  particular 
filter  in  the  look  up  table  since  one  filter  may  be 
suitable  for  different  ship  characteristics  and  sea 
conditions . 

•  The  weighting  factor  A  used  in  the  performance 
criterion  equation  3.5  plays  an  exceptional  and  impor- 
tant role  in  the  optimal  controller  parameters  determi- 
nation.  However,   it  is  obtained  from  studies  based  on 
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many  assumptions  and  therefore  isn't  predicted  accu- 
rately. As  a  consequence  the  controller  found  does  not 
minimize  added  resistance  unless  the  proper  value  of 
the  weighting  factor  A  has  been  used. 

•  The  method  used  to  obtain  the  average  for  the  added 
mass  and  added  inertia  for  the  irregular  seas  studies 
might  not  represent  the  actual  average. 

•  The  function  minimization  subroutine  BOXPLX  after  two 
modifications  seems  to  be  pretty  suitable  for  on  board 
use  working  as  a  main  part  of  the  adaptive  scheme. 

B.   RECOMMENDATIONS  FOR  FUTURE  STUDIES 

The  following  recommendations  for  future  work  to  gain  a 
fuller  and  deeper  understanding  of  the  problem  are  made  as 
follows : 

•  Some  studies  are  necessary  in  order  to  investigate  why 
part  of  the  obtained  controllers  in  Tables  V,  VI,  VII, 
VIIT  are  lag  and  part  of  them  are  lead  filters. 

•  As  we  stated  in  chapter  4  the  yaw  and  rudder  excursions 
in  Figures  4.2  through  4.9  are  less  than  1*.  It  is 
necessary  to  investigate  the  reason  for  that.  It  might 
be  because  of  the  optimization  of  the  filter  or  the 
forces  and  moments  have  not  been  sealed  properly. 

•  It  is  necessary  to  find  out  the  appropriate  average 
values  for  the  added  mass  and  added  inertia  before  we 
attempt  further  studies  in  irregular  seas. 

•  The  full  hydrodynamic  coefficients  for  the  SL-7  are 
necessary  in  order  to  develop  and  include  the  surge 
equation  in  our  ship  model.  So  far  we  ignore  the  surge 
equation  and  we  assumed  constant  ship  speed  while  in 
reality  the  added  resistance  due  to  steering  must 
reduce  the  ship  speed. 
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•  By  developing  the  surge  equation  in  our  ship  model  we 
may  be  able  to  determine  a  good  value  for  the  weighting 
factor  to  be  used  in  the  performance  criterion. 

•  Also,  with  the  surge  equation  available  in  our  model  we 
should  be  able  to  calculate  actual  energy  losses  and 
savings  in  fuel. 

•  As  we  stated  in  chapter  6  the  time  frames  (samples) 
used  by  the  BOXPLX  must  be  long  relative  to  the  time 
required  for  the  initial  condition  response  to  die  out. 
Reid  [Ref.  1,2]  has  chosen  time  frame  600  seconds  long. 
This  time  frame  is  long  and  it  is  necessary  to  find  out 
the  sufficient  time  frame  since  the  time  required  by 
the  function  minimization  subroutine  for  convergence  is 
directly  proportional  to  the  sample  duration. 
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APPENDIX  A 
NOMOTO  THIRD  ORDER  MODEL  DETERMINATION 

//NOMOTO  JOB  (XXXX,X.XXX)  ,  'RESEARCH'  ,  CLASS  =  J 

//-MAIN  ORG=NPGVMl.XXXXP 

//  EXEC  FORTXCG,PARM.FORT= 'OPT (2) ' , IMSL=DP ,REGION= 1024K 

//FORT. SYS IN  DD  * 

C   IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 

C   OBTAINED  CHANGE  XS(*)  TO  X(*)  AND  DELETE  XU  (•'•),  AND  XL  ( -  ) 

DIMENSION  XS(4) ,XU(4) ,XL(4) 

XS(1)=0. 1 

XS(2)=15. 13 

XS(3)=15.675 

XS(4)=9.014 
C   XS(I)  IS  THE  STARTING  GUESS 

C   XL(I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 
C  •  XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I ' TH  VARIABLE 

XL(1)=.01 

XU ( 1 )  =  1 . 0 

XL(2)=1.0 

XU(2)=20. 

XL ( 3 ) = L . 0 

XU(3)=100. 

XL(4)=1.0 

XU(4)=100.0 
C   A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C   IS  DISCUSSED  IN  BOXPLX 

R=9./13. 

NTA=1000 

NPR=0 

NAV=0 

NV  =  4 

IP=0 
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C   THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 

C   CALL  PLANT (XX) 

C   IF  ONLY  SIMULATION  IS  WANTED 

CALL  BOXPLX ( NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 

WRITE  (6,25) 
25    FORMAT (IX,'  OPTIMAL  GAINS',/) 

DO  30  1=1,4 
30    WRITE(6,40)I,XS(I) 
40    F0RMAT(1X, 'X(' ,12, ' )=' ,F14.7) 

WRITE  (6,77)  TDIFF 
77    F0RMAT(1X, ' COST= ' ,F14.7) 

STOP 

END 
C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
-  SUBROUTINE  PLANT (XX) 

COMMON  TDIFF 

REAL*8  L,L2,L3 ,L4,L5 ,L6 , RXR , RXI , RYR ,RYI ,MZR ,MZI , RX , RY 

REAL -8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT , TX , TY 

REAL*8  TIME , ETIME , XUDOT , XU , XUU , XVR , XVV , XDD , WA , WE , RZ 

REAL- 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT , TZ 

REAL- 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL- 8  RHO , IZ , FX , FY , MZ , XP , MAS S , DELT 

REAL- 8  YAWE,YAWC,D2,D 

REAL- 8  K , TP 1 , TP2 , Z , XI , X2 , X 3 , X4 , X5 , YAW2 , S 

DIMENSION  XX (4) 
C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 

ETIME=600. 

TIME=0.0 

ICOUNT^l 
C   INITIALIZE  THE  COST  FUNCTION 

TDIFF=0.0 
C   GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 

K=XX(1) 

Z=XX(2) 
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TP1=XX(3) 

TP2=XX(4) 
C      X,XDOT,Y,YDOT    ARE    FIX    COORDINATES    ON    EARTH 

X  =  0.0 

Y=0.0 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW=0.0 

R=0.0 

RDOT=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C   38.81  FT/SEC=23  KNOTS 

UC=38.81 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z  ' 

FX=0.0 

FY=0.0 

MZ=0.0 
C      RXR=- . 15744D+05 
C      RXI=- . 19950D+06 
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C  RYR=0.52365D+04 
C  RYI=0. 18699D+06 
C  MZR=- .29870D+08 
C      MZI=- .37751D+07 

RXR=- .50230D+04 

RXI=0.12712D+05 

RYR=0.35290D+04 

RYI=- .31909D+05 

MZR=0.38826D+07 

MZI=- .64313D+07 
C  RXR=0.28540D+04 
C  RXI=- .99574D+04 
C  RYR=- .85441D+04 
C  RYI=0.39595D+05 
C  MZR=-.13014D+08 
C  MZI=0.11348D+08 
C  RXR=- .75642D+04 
C  RXI=0.83497D+04 
C  RYR=0.23379D+05 
C  RYI-- .81502D+05 
C  MZR=0.28622D+07 
C  MZI=-.19388D+08 
C  RXR=-.37916D+04 
C  RXI=0. 16381D+04 
C  RYR=- . 76647D+05 
C  RYI=- .37685D+05 
C  MZR=- .83915D+07 
C      MZI=-.53176D+07 

RX=DSQRT (RXR**2+RXI**2 ) 

RY  =  D  S QRT ( RYR - * 2  +  R YI ** 2 ) 

RZ=DSQRT (MZR**2+MZI**2 ) 

TX = DATAN2 ( RXI , RXR ) 

TY=DATAN2 ( RYI , RYR ) 

TZ=DATAN2 (MZI ,MZR) 
C   SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-0.32,2-0.75,3-2.5, 
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C   4-5.0,5-7.0,6-10.0,7-17.5,8-20.5,9-27.0 

WA=10.0 
C   ENCOUNTER  FREQUENCY  .1,. 2,  .3,  .4,  .5,  .6,  .75, 1.0, 1.5, 2. 5 

WE=0.4 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 

RHO=1.9876 

MASS= ( . 0044) *( . 5*RHO*L3 ) 

IZ= (0 . 00028)*( . 5*RHO*L5 ) 

YAWE=0.0 

X1=0.0 

X2=0.0 

X3=0.0 

X4=0.0 

X5=0.0 

YAW2=0.0 
200  CONTINUE 

S  =  DSQRT(U"-'2  +  V""2) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC= ( 1 . 0/ 57 . 296 ) 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER  (YAW  ACTUAL  -  YAW  COMMAND) 
C   FOR  EQUATIONS  OF  MOTION. 

YAWE=YAW  -  YAWC 

D=YAWE 
C 

C     NOMOTO  3RD  ORDER  PLANT 
C 

C   ERROR  SIGNAL  TO  DRIVE  RUDDER  (YAW  COMMAND  -  YAW  ACTUAL) 
C   FOR  NOMOTO  MODEL. 

D2=YAWC-YAW2 

X1=(D2-X2)/TP1 

X3=K*(Z*X1+X2) 

X4=(X3-X5)/TP2 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

XUDOT= ( - . 000 1 ) * ( . 5 *RHO*L3 ) 

69 


XUU= ( -  0 . 0003 ) * ( . 5 -'RHO-L2 ) 

XU= ( -  0 . 025  3 ) * ( 0 . 5 -RHO-L2 * S ) 

XVR=  ( 0  .  00 3  9  )«'-(.  5  -RHO "L3  ) 

XVV  = ( - . 0012 )*( . 5 "RHO- L 2 ) 

XDD  = ( - 0 . 0005 ) * ( . 5 -RHO *L2 * S * * 2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
C      YV= ( -  0 . 00  7  5  8 ) * ( . 5 -RHO -L2" S ) 

YR=  ( 0  .  0023  )  -  (' .  5 -RHO -L3  -  S  ) 

YD  = ( 0 . 00 145 ) - ( . 5 -RHO -L2 - S - -  2 ) 

YVVR= ( 0 . 0 1 ) - ( . 5 -RHO - L3 / S ) 

YVRR= ( -  0 . 008 ) * ( . 5 -RH0-L4/ S ) 

YVVV= ( -  0 . 0  3 ) - ( • 5  - RHO - L2 / S ) 

YRRR= ( 0 . 00  3 ) * ( . 5 -RHO -L5 / S ) 

YDDD= ( -0 . 0005 )- ( . 5 -RHO-L2- S * -2 ) 

YVDOT=-0.30908D+07 

YV=-0.81271D+04 
C      YVDOT=- .36185D+07 
C      YV=- .24757D+06 
C      YVDOT=- .  32890D+07 
C      YV=- .11775D+07 
C      YVDOT=- .23038D+07 
C     YV=-.18267D+07 
C     YVDOT=- .59800D+06 
C     YV=-.13260D+07 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= ( -  0 . 00213 )*( . 5 -RHO-L3 * S ) 
C      NR= (-0.00105) * ( . 5 -RH0-L4-S ) 

ND= ( -0 . 0007 )*( . 5-RH0-L3-S- -2 ) 

NVVR= ( -  0 . 0 15 ) - ( • 5 -RHO - L4 / S ) 

NVRR= ( -  0 . 0  0  8 ) - ( . 5  * RHO -L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5 -RHO -L3 / S ) 

NRRR= ( -  0 . 006 ) * ( . 5 -RH0-L6 / S ) 

NDDD  = (0.0001) - ( . 5 -RHO -L3 - S - -2 ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
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c 

C     NRDOT  =  ( - 0  .  00027  )'*  (  .  5 *RHO *L5  ) 

C   SPEED=23  KNOTS, ENCOUNTER  ANGLE  =  60, ENCOUNTER  FREQ=0.75 
C      NRDOT=- .26251D+12 
C      NR=- .53637D+09 

NRDOT=- .20125D+12 

NR=-  .94970D+10 
.C      NRDOT=- . 18671D+12 
C      NR=- .46860D+11 
C      NRDOT=-. 14518D+12 
C      NR=- .87538D+11 
C      NRDOT=- .37261D+11 
C      NR=- .69856D+11 
C   SETS  SEA  STATE  TO  ZERO 

FX=0. 

FY  =  0. 

MZ  =  0. 

FX=WA*RX*DCOS (WE*TIMF>TX ) 

FY  =  WA* RY *DCO S ( WE * TIME  +  TY ) 

MZ=WA*RZ*DCOS (WE-TIME+TZ ) 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU"UC""2 
C   EQUATIONS  OF  MOTION 

C      UDOT=(  (XVR  +  MASS)*V*R  +  XUU*U**2  +  XVV*V**2 
C     1  +  XDD-D-D  +   FX  +  XP  ) / (MASS-XUDQT ) 

VDOT=(YV*V  +  (YR-MASS*U)*R  +  YD*D  +  YWR*V**2*R 

1  +  YVRR*V*R**2  +  YVW*V*V**3 

2  +  YRRR*R**3  +  YDDD*D**3  +  FY  ) / (MASS- YVDOT ) 
RDOT=(NV*V  +  NR-R  +  ND*D  +  NVVR*V**2*R  +  NVRR*V*R**2 

1  +NWV*V**3  +  NRRR*R**3  +  NDDD-D--3  +  MZ  )/(IZ- NRDOT) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 
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C   CONVERT  RADIANS  TO  DEGREES 
5  0   YAWDEG  =  YAW- 5  7.296 

RDEG=R-5  7.29  6 

RDDEG=RDOT-5  7 .29  6 

DDEG=D-5  7.29  6 

YAWC=YAWC-5  7.29  6 

ICOUNT=l 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1.0 
C   INTEGRATION 

X2=X2+X1-DELT 

X5=X5+X4-DELT 

YAW2=YAW2+X5-DELT 

U=U+UDOT*DELT 

V=V+VDOT*DELT 

R=R+RDOT*DELT 

YAW=YAW+R*DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 

XDOT  =  U * DCO S  ( YAW )  -  V*D-SIN  ( YAW ) 

YDOT=U*DSIN(YAW)+V*DCOS(YAW) 

X=X+XDOT*DELT 

Y=Y+YDOT*DELT 

TIME=TIME+DELT 

ICOUNT=ICOUNT+l 
C   COST  FUNCTION 

TDIFF=TDIFF*  (YAW-YAW2  ) —2 

GO  TO  200 
400   CONTINUE 
C      WRITE(6,500)  TDIFF ,K , Z , TP1 , TP2 
500   FORMAT ('  ' , IX , '   COST  = ' ,F12 . 7 , 2X , '  K  =',F10.7, 
I  '  Z  =',F15.7S'  TP1  =' ,FL5.7,2X, '  TP2  =',F15.7) 

RETURN 

END 

72 


APPENDIX  B 
REGULAR  SEASTATE  FORMULATION 

//REGUGAINS  JOB  (XXXX , XXXX) ,' RESEARCH ', CLASS=C 
//-MAIN  ORG=NPGVMl.XXXXP 

//  EXEC  FORTXCG,PARM.FORT= TOPT(2) ' , IMSL^DP , REGION^ 1024K 
//FORT. SYS IN  DD  * 

C   IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 
C   OBTAINED  CHANGE  XS(*)  TO  X(*)  AND  DELETE  XU(*) ,AND  XL(-) 
DIMENSION  XS(3) ,XU(3) ,XL(3) 
XS(1)=0. 9650610 
XS(2)=0. 4500911 
XS(3)=5. 6194260 
C   XS(I)  IS  THE  STARTING  GUESS 

C   XL(I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 
C   XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I'TH  VARIABLE 
.  XL(1)=0.1 
XU(1)=4.0 
XL(2)=0.1 
XU(2)=15.0 
XL(3)=1.0 
XU(3)=25.0 
C   A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C   IS  DISCUSSED  IN  BOXPLX 
R=9./13. 
NTA=1000 
NPR=100 
NAV=0 
NV=3 
IP=0 
C   THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C   CALL  PLANT (X) 
C   IF  ONLY  SIMULATION  IS  WANTED 
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CALL  BOXPLX (NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
WRITE  (6,25) 

25      FORMAT (IX,'  OPTIMAL  GAINS',/) 
DO  30  1=1,3 

30      WRITE(6,40)I,XS(I) 

40      F0RMAT(1X,  'X( '  ,12,  ' )  = '  ,F14.7) 
STOP 
END 
SUBROUTINE  PLANT (XX) 

C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL- 8  L,L2,L3 ,L4,L5 ,L6 

REAL -8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL- 8  TIME , ETIME , XUDOT , XUU , XVR , XW , XDD 
REAL  -8  YV  ,  YR  ,  YD  ,  YWR ,  YVRR  ,  YVW  ,  YRRR ,  YDDD  ,  YVDOT 
REAL -8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL -8  RHO , IZ , FX , FY , MZ , XP , MAS S , DELT , MZI , WA , WE 
REAL- 8  DYAWE , YAWE , YAWC , ISE , ISR , LAMDA , D , RYR , RYI , MZR 
REAL- 8  Kl , Tl , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR , RXI 
DIMENSION  XX(3) 

C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 

C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 

C   SIMULATION  END  TIME  IN  SECONDS 
ETIME=600.0 
TIME=0.0 
ICOUNT=l 

C   INITIALIZE  THE  COST  FUNCTION 
ISE=0.0 
ISR=0.0 
TDIFF=0.0 
LAMDA=8 .128 

C-  GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K1=XX(1) 
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T1=XX(2) 
T2=XX(3) 

C      WRITE(6,1010)  K1,T1,T2 

C1010   FORMAT (IX, *K1  =',F15.7,'  Tl  =',F15.7,'  T2  =',F15.7) 

C   X,XDOT,Y,YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X=0.0 

Y=0.0 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT^O.O 

VDOT=0.0 

YAW=0.0 

R=0.0 

RDOT^O.O 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.81  FT/SEC=23  KNOTS 

UC=38 .81 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 

75 


C  RXR=-0.91037D+03 

C      RXI=0.50869D+05 

C      RYR=-0.20256D+04 

C      RYI=0. 18077D+06 

C      MZR=- . 14310D+08 

C      MZI=- . 16903D+07 

C     RXR=-0.99047D+04 

C      RXI= . 15994D+06 

C      RYR=- .64455D+05 

C      RYI=0.61873D+06 

C      MZR= . 120180+08 

C      MZI=- .49204D+07 

C  RXR=-0.32876D+05 

C      RXI=0.25844D+06 

C      RYR=-.27053D+06 

C      RYI=.90191D+06 

C      MZR=0. 11964D+09 

C     MZI=0.24103D+08 

C      RXR=- .54639D+05 

C      RXI=.28236D+06 

C      RYR=- .28668D+06 

C      RYI=0.79670D+06 

C     'MZR=0.19925D+Q9 

C      MZI=0.77746D+08 
RXR=0.27268D+05 
RXI=- .71601D+05 
RYR=0.14077D+05 
RYI  =  - .28679D+06 
MZR=- . 30892D+08 
MZI=- .53246D+08 
RX  =  D  S QRT ( RXR * * 2  +  RX I * * 2 ) 
RY  =  D  S QRT ( RYR* * 2  +  R Y I * * 2 ) 
RZ=DSQRT (MZR**2+MZI**2 ) 
TX = D ATAN2  ( RX I , RXR ) 
TY= DATAN2 ( RYI , RYR ) 
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TZ=DATAN2(MZI,MZR) 
C   SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-0.32,2-0.75,3-2.5, 
C   4-5. 0,5-7. 0,6- 10. 0,7-17. 5, 8-20. 5, 9-27.0 

WA=10.0 
C   ENCOUNTER  FREQUENCY  .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 

WE=1.5 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 

RHO=1.9876 

MASS  = ( . 0044 ) * ( . 5 -RHO-L3 ) 
-IZ=(0.00028)*(.5 * RHO * L5 ) 

YAWE=0.0 

X2=0.0 

DX2=0.0 
200  CONTINUE 

S=DSQRT(U**2+V**2 ) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C.  ERROR  SIGNAL  TO  DRIVE  RUDDER(YAW  ACTUAL  -  YAW  ORDERED) 
C   (  COMPENSATOR  FILTER  ) 

YAWE=YAW  -  YAWC 

DX2=(YAWE-X2)/T2 

D=K1*(T1*DX2+X2) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

XUDOT = ( - . 000 1 ) * ( . 5 -RHO -L3 ) 

XU= ( -0 . 0253 )*( . 5*RHO*L2*S ) 

XUU  =(-0.0003) * ( . 5 * RHO * L2 ) 

XVR= ( 0 . 0039 )* ( . 5-RHO-L3 ) 

XVV= ( - . 00 12 ) * ( . 5 -RHO-L2 ) 

XDD  = ( -  0 . 0  0  0  5 ) * ( . 5 -RHO - L2  * S  ** 2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
C      YV= ( -  0 . 00  7  5  8 ) - ( . 5 -RHO -L2 * S ) 
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YR  = ( 0 . 0  0  2  3 ) * ( . 5 * RHO  * L  3 * S ) 

YD= ( 0 . 00 145 ) * ( . 5 -RHO-L2 - S * * 2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 -RHO -L3 / S ) 

YVRR= ( - 0 . 008 ) * ( . 5 * RHO * L4 / S ) 

YVVV= ( - 0 . 03 ) * ( . 5 * RHO *L2 / S ) 

YRRR= ( 0 . 003 ) * ( . 5 * RHO *L5 / S ) 

YDDD= ( -  0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

C      YVDOT= ( -  0 . 003  9 ) * ( . 5*RHO*L3 ) 

C   SPEED=23  KNOTS, ENCOUNTER  ANGLE= 60 , ENCOUNTER  FREQ=0.75 
C      YVDOT=- . 30908D+07 
C      YV=- .81271D+04 
C      YVDOT=-.36185D+07 
C      YV=-.24757D+06 
C      YVDOT=- .32890D+07 
C      YV=- .11775D+07 
C      YVDOT=- .23038D+07 
C      YV=- . 18267D+07 

YVDOT=- .59800D+06 

YV=- .13260D+07 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 


NV=(-0. 00213)* 
NR= (-0.00 105) * 
ND=(-0.0007)*( 
NWR=(-0.015)* 
NVRR= (-0.008 )* 


NVVV  = ( 0 . 0 1 ) * ( . 5 * RHO *L3 / S ) 


NRRR= (-0.006)* 
NDDD= (0.0001)* 


.5*RHO*L3*S) 

.5-RHO-L4-S) 
5"RHO"L3"S""2) 
.5*RHO*L4/S) 
.5-RHO-L5/S) 


5-RHO-L6/S) 
5*RHO*L3*S**2) 

C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 

C   FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 

C 

C      NRDOT  =  (  -  0  .  0  0  0  2  7  )  *  (  .  5  *  RHO  *  L5  ) 
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C   SPEED=23  KNOTS, ENCOUNTER  ANGLE= 60 , ENCOUNTER  FREQ=0.75 

C      NRDOT=- .26251D-12 

C      NR=- .53637D+09 

C      NRDOT=- .20125D+12 

C      NR=- . 94970D+10 

C      NRDOT=- . 18671D+12 

C      NR=- .46860D+11 

C      NRDOT=- . 14518D+12 

C      NR=- .87538D+11 

NRDOT=- .37261D+11 

NR=- .69856D+11 
C   REGULAR  WAVE  SEA  STATE 

FX  =  WA*RX*DCO  S ( WE  - T IME  +  TX ) 

FY  =  WA-RY-D  CO  S ( WE - T IME  +  TY ) 

MZ=WA*RZ*DCOS (WE*TIME+TZ ) 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU*UC**2 
C   EQUATIONS  OF  MOTION 

C      UDOT=(  (XVR  +  MASS)*V*R  +  XUU*U**2,  +  XW*V**2 
C     1  +  XDD*D*D  +   FX  +  XP  ) / (MASS-XUDOT) 

VDOT=(YV-V  +  (YR-MASS-U)-R  +  YD*D  +  YWR*V**2*R 

1  +  YVRR*V*R**2  +  YVW*V**3 

2  +  YRRR*R**3  +  YDDD*D**3  +  FY  ) / (MASS- YVDOT) 
RDOT=(NV*V  +  NR-R  +  ND*D  +  NWR*V**2 *R  +  NVRR*V*R**2 

1  +  NVW*V**3  +  NRRR*R**3  *  NDDD-D--3  +  MZ )  /  (IZ-NRDOT) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
50   YAWDEG=  YAW -5 7. 29 6 
RDEG=R*57.296 
RDDEG=RDOT*57.296 
DDEG=D-5  7.29  6 
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YAWC=YAWC*5  7 .29  6 

IC0UNT=1 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1.0 
C   INTEGRATION 

U=U+UDOT-DELT 

V=V+VDOT*DELT 

R=R+RDOT*DELT 

YAW=YAW+R*DELT 

X2=X2+DX2*DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C      XDOT=U*DCOS  ( YAW )  - V-DSIN ( YAW ) 
C      YDOT=U*DSIN(YAW)+V*DCOS(YAW) 
C      X=X+XDOT*DELT 
C      Y=Y+YDOT*DELT 

TIME = TIME + DELT 

I COUNT = I COUNT +1 

ISE=ISE  +  LAMDA*YAWE**2 

ISR=ISR  +  D**2 

GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE (6, 5 00)  ISE,ISR,TDIFF,K1,T1,T2 
500   FORMAT ('  ' , IX , ' ISE= ' , F15 . 7 , '   ISR= ' , F15 . 7 , '   TOTAL= ' , 
1F15.7,2X,  'Kl=  '  ,F15.7,2X,  *T1=  '  ,F15  .7  ,  2X/T2  =  '  ,F15.7) 

RETURN 

END 
The  function  minimization  subroutine  BOXPLX  follows. 
Then  the  following  two  cards  must  be  placed. 
//GO.SYSIN  DD  * 
/* 
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APPENDIX  C 
SYSTEM'S  RESPONSE  FOR  REGULAR  SEAS 

//REGURESP  JOB  (XXXX ,XXXX) ,  ' RESEARCH ', CLASS=A 
//-MAIN  ORG=NPGVMl.XXXXP 

//  EXEC  FORTXCG,PARM.FORT= 'OPT(2) ' ,IMSL=DP,REGION=1024K 
//FORT.SYSIN  DD  * 

C   IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 
C   OBTAINED  CHANGE  XS(")  TO  X(*)  AND  DELETE  XU(-),AND  XL(*) 
COMMON  J 
DIMENSION  X(3) 
X(l)=l. 5420017 
X(2)=141. 2350922 
X(3)=23. 8943634 
C   CALL  PLANT (X) 

C   IF  ONLY  SIMULATION  IS  WANTED 
CALL  PLANT (X) 
WRITE  (6,25) 
25     FORMAT (IX,'  OPTIMAL  GAINS',/) 

DO  30  1=1,3 
30     WRITE(6 ,40)I,X(I) 
40     F0RMAT(1X, 'X( ' ,12, ' )= ' ,F14.7) 
WRITE (6, 50)  J 
50     F0RMAT(1X,'J  =  ',E15.10) 
STOP 
END 

SUBROUTINE  PLANT (XX) 
C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL- 8  L,L2,L3 ,L4,L5 ,L6 

REAL- 8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL- 8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL- 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
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REAL- 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL -8  RHO , I Z , FX , FY , MZ , XP , MAS  S , DELT , MZI , WA , WE 

REAL- 8  DYAWE , YAWE , YAWC , ISE , ISR , LAMDA , D , RYR , RYI , MZR 

REAL-8  Kl , Tl , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR , RXI 

DIMENSION  XX(3) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 

ETIME=600. 

TIME^O.O 

ICOUNT=l 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=8. 128 
C   GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 

K1=XX(1) 

T1=XX(2) 

T2=XX(3) 
C      WRITE(6,1010)  K1,T1,T2 

C1010   F0RMAT(1X, 'Kl  =',F15.7,'  Tl  =f,F15.7,'  T2  =',F15.7) 
C   X,XDOT,Y,YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X  =  0.0 

Y=0.0 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW=0.0 

R=0.0 
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RDOT=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.82  FT/SEC=23  KNOTS 

UC=38.82 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 

RXR=- . 91037D+03 

RXI=0.50896D+05 

RYR=- .20256D+04 

RYI  = . 18077D+06 

MZR=- . 14310D+08 

MZI=- . 16903D+07 

RX  =  DSQRT ( RXR* -  2  +  RXI * - 2 ) 

RY  =  D  S QRT ( RYR - - 2  +  RY I * * 2 ) 

RZ  =  D  S QRT ( MZR - - 2  + MZ I ** 2 ) 

TX=DATAN2 ( RXI , RXR ) 

TY=  DATAN2 ( RYI , RYR ) 

TZ=DATAN2(MZI,MZR) 
C   SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-0.32,2-0.75,3-2.5, 
•C   4-5.0,5-7.5,6-10.0,7-17.5,8-20.5,9-27.0 

WA=27.0 
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C   ENCOUNTER  FREQUENCY  .1,. 2,  .3,  .4,  .5,  .6,  .75, 1.0, 1.5, 2. 5 

WE=0.2 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 

RHO=1.9876 

MAS  S  = ( . 0044 ) * ( . 5 -RHO -L3 ) 

IZ=(0.00028)*(-5 *RHO - L5 ) 

YAWE=0.0 

X2=0.0 

DX2=0.0 
200  CONTINUE 

S  =  DSQRT (U* *2  +  V- - 2 ) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 
C   (  COMPENSATOR  FILTER  ) 

YAWE=YAW  -  YAWC 

DX2=(YAWE-X2)/T2 

D=K1"(T1-DX2+X2) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

XUDOT= ( - . 0001 ) * ( . 5 -RHO-L3 ) 

XU= ( - 0 . 025  3 ) * ( . 5 *RHO*L2*S ) 

XUU=  (  -  0  .  0 0 0  3  )  *  (  .  5  * RHO * L2  ) 

XVR= (0.0039) * ( . 5 -RHO *L3 ) 

XW=  (  -  .  00 12  )  *  (  .  5*RHO*L2  )  ' 

XDD= ( -  0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
C      YV=(-0.00758)-( .5-RHO-L2-S) 

YR  = ( 0 . 0  0  2  3 ) * ( . 5 -RHO  * L  3 * S ) 

YD= ( 0 . 00 145 ) * ( . 5 -RHO-L2 - S -  - 2 ) 

YVVR= ( 0 . 0 1 ) * ( , 5 -RHO * L3 / S ) 

YVRR= ( -  0 . 0 0 8 ) - ( • 5 -RHO - L4 / S ) 
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YVVV= ( -  0 . 03 ) * ( . 5 -RHO -L2 / S ) 

YRRR= ( 0 . 0  0  3 ) * ( . 5 * RHO * L5 / S ) 

YDDD  = ( -  0 . 0005 ) * ( . 5 *RHO*L2 *S * *2 ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

C      YVDOT  = ( - 0 . 00  3  9 ) * ( . 5 * RHO * L3 ) 
C   SPEED=23  KNOTS, ENCOUNTER  ANGLE= 60 , ENCOUNTER  FREQ=0.75 

YVDOT=- .30908+07 

YV=-0.81271D+04 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= ( -  0 . 00213 ) * ( . 5*RHO*L3*S ) 
C      NR= ( -  0 . 0  0 1 0  5 ) * ( . 5 -RHO  *L4* S ) 

ND  = ( -  0 . 0 0 0 7 ) * ( . 5 -RHO * L 3 * S ** 2 ) 

NVVR= ( -0 . 015 ) * ( . 5 -RH0-L4/ S ) 

NVRR= ( -  0 . 0  0 8 ) * ( . 5 * RHO * L5 / S ) 

NVVV= ( 0 . 0 1 ) - ( . 5 -RHO -L3 / S ) 

NRRR= ( -  0 . 0 0 6 ) - ( . 5 -RHO * L 6  /  S ) 

NDDD  = ( 0 . 0  0  0 1 ) - ( . 5  * RHO  * L  3  - S - -  2 ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
C 

C      NRDOT  = ( -  0 . 0 0 0 2  7 ) - ( . 5  * RHO * L5  ) 
C   SPEEDS  32  KNOTS, ENCOUNTER  ANGLE= 60 , ENCOUNTER  FREQ=0.75 

NRDOT=-0.26251D+12 

NR=-0.53637D+09 
C   REGULAR  WAVE  SEA  STATE 

FX  =  WA"RX-DCOS (WE-TIME  +  TX) ' 

FY=WA"RY-DCOS (WE-TIME+TY ) 

MZ=WA-RZ"DCOS (WE-TIME+TZ ) 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU-UC""2 
C   EQUATIONS  OF  MOTION 
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C      UDOT=(  (XVR  +  MASS)*V*R  +  XUU*U**2  +  XW*V**2 
C     1  +  XDD-D-'D  +   FX  +  XP  )  /  (MASS-XUDOT ) 

VDOT=(YV*V  +  (YR-MASS-U)---R  +  YD*D  +  YVVR*V**2*R 

1  +  YVRR-V-R--2  +  YVW*V**3 

2  +  YRRR-R— 3  +  YDDD-D--3  +  FY  ) / (MASS -YVDOT ) 
RDOT=(NV*V  +  NR*R  +  ND*D  +  NVVR*V**2*R  +  NVRR*V*R**2 

1  +  NVW*V**3  +  NRRR*R**3  +  NDDD*D**3  +  MZ )  /  ( IZ-NRDOT ) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.  2)  GO  TO  50 

GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
5  0   YAWDEG  =  YAW -5  7.296 

RDEG=R*57.296 

RDDEG=RDOT*57.296 

DDEG=D*57.296 

YAWC=YAWC*5  7 .29  6 

WRITE  (6,101)  TIME, YAWDEG 
101   FORMAT ( IX, F 12. 8 ,1X,F12.8) 

I COUNT =1 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1.0 
C   INTEGRATION 

U=U+UDOT*DELT 

V=V+VDOT-DELT 

R=R+RDOT*DELT 

YAW=YAW+R*DELT 

X2=X2+DX2*DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C      XDOT=U-DCOS(YAW)-V"DSIN(YAW) 
C      YDOT=U*DSIN ( YAW ) + V-DCOS ( YAW ) 
C      X=X+XDOT-DELT 
C      Y=Y+YDOT*DELT 

TIME=TIME+DELT 
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IC0UNT=IC0UNT+1 
ISE^ISE  +  LAMDA*YAWE**2 
ISR=ISR  +  D**2 
GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE (6, 5 00)  ISE,ISR,TDIFF,K1,T1,T2 
500   FORMAT C  '  , IX ,  ' ISE= '  , F15 . 7 ,  '   ISR= '  , F15 . 7  ,  '   TOTAL= T  , 
1  F15.7,2X, 'Kl=' ,F15.7,2X, 'Tl=' ,F15.7,2X, 'T2= ' ,F15.7) 
RETURN 
END 
//GO. SYS IN  DD  * 

/* 
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APPENDIX  D 
IRREGULAR  SEASTATE  FORMULATION 

//IRREGAINS  JOB  (XXXX ,XXXX) ,' RESEARCH ', CLASSIC 

//-MAIN  ORG=NPGVMl.XXXXP 

//  EXEC  FORTXCG,PARM.FORT='OPT-(2)  '  ,IMSL  =  DP,REGION=1024K 

//FORT. SYS IN  DD  * 

C   IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 

C   OBTAINED  CHANGE  XS(*)  TO  X(*)  AND  DELETE  XU(*),AND  XL(-) 

DIMENSION  XS(3) , XU ( 3 ) , XL ( 3 ) 

XS(1)=0. 655751 

XS(2)=90.5483 

XS(3)=36. 74847 
C   XS(I)  IS  THE  STARTING  GUESS 

C   XL (I)  IS  THE  LOWER  LIMIT  FOR  THE  I  * TH  VARIABLE 
C   XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I ' TH  VARIABLE 

XL  ( 1 )  =  0  . 0 1 

XU ( 1 )  =  2  .  0 

XL  (.2)  =  20.0 

XU(2)=180.0 

XL ( 3  )  =  5  .  0 

XU(3)=180.0 
C   A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C   IS  DISCUSSED  IN  BOXPLX 
R=9./13. 
NTA=1000 
NPR=100 
NAV=0 
NV=3 

ip=6 

C   THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 

C   CALL  PLANT (X) 

C   IF  ONLY  SIMULATION  IS  WANTED 
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CALL  BOXPLX (NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
WRITE  (6,25) 
25      FORMAT (IX,'  OPTIMAL  GAINS ',/ ) 

DO  30  1=1,3 
30      WRITE(6,40)I,XS(I) 
40      F0RMAT(1X, 'X( ' ,12, ' )= ' ,F14.7) 
STOP 
END 

SUBROUTINE  PLANT (XX) 
C   SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 

COMMON  TDIFF 

REAL-8  L,L2,L3 ,L4,L5 ,L6 

REAL -8 . X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 

REAL- 8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 

REAL- 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 

REAL- 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL- 8  RHO , IZ , FX , FY , MZ , XP , MASS , DELT 

REAL -8  DYAWE , YAWE , YAWC , ISE , ISR , LAMDA , D 

REAL-8  K1,TI,T2,T3 , T4 , D , X2 , DX2 , X3 ,DX3 ,X4,CH(11) ,S 

DIMENSION  XX(3) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

C   INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 

ETIME=600. 

TIME=0.0 

ICOUNT=l 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=8 . 128 
C   GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 

K1=XX(1) 
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T1=XX(2) 

T2=XX(3) 
C  WRITE(6,1010)    K1,T1,T2 

C      X,XDOT,Y,YDOT    ARE    FIX    COORDINATES    ON    EARTH 

X=0.0 

Y=0.0 

XDOT=0.0 

YDOT=0.0- 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW=0.0 

R=0.0 

RDOT=0.0 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.82  FT/SEC=23  KNOTS 

UC=38.82 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2  =  L-"V2 

L3=L*L-'L 

L4=L*L3 

L5=L*L4 

L6=L-L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 
C   ISEA  IS  A  SWITCH  ;ISEA=0  (CALM  WATER)  ISEA=1  (SEA  STATE) 

90 


ISEA=1 

C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 

RHO=1.9876 

MASS= ( . 0044 ) * ( . 5*RHO*L3 ) 

IZ= ( 0 . 00028 ) * ( . 5 *RHO *L5 ) 

YAWE=0.0 

X2=0.0 

DX2=0.0 
200  CONTINUE 

S  =  D  S  QRT  ( U  *  *  2  +  V  *  *  2  ) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  ( TIME. GE. 0.0)  YAWC=0.0 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 
C   (  CONTROLLER   FILTER  ) 

YAWE=YAW  -  YAWC 

DX2=(YAWE-X2)/T2 

D=K1*(T1*DX2+X2) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
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XUDOT = ( - . 000 1 ) * ( . 5*RHO*L3 ) 
XU= ( - 0 . 025  3 ) * ( . 5 *RHO*L2*S ) 

XUU= ( -  0 . 000  3 ) * ( . 5 -RHO-L2 ) 
XVR= ( 0 . 0  0  3  9 ) * ( . 5 * RHO * L3 ) 
XVV= ( - . 0012 ) * ( . 5*RHO*L2 ) 
XDD=  (  -  0  .  0005  )  *  (  .  5  *RHO*L2  *  S  *  *  2  ) 
LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
"  YV= ( -  0 . 00758 ) * ( . 5*RHO*L2*S ) 
YR= ( 0 . 0023 )* ( . 5 *RHO*L3 *S ) 
YD= ( 0 . 00 145 ) - ( . 5 -RHO *L2 *S - * 2 ) 
YVVR= ( 0 . 0 1 ) * ( . 5 * RHO »L3 / S ) 
YVRR= ( -  0 . 0  0  8 ) * ( . 5 * RHO * L4 / S ) 
YVVV  = ( - 0 . 0  3 ) * ( . 5  - RHO - L2 / S ) 
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YRRR  = ( 0 . 003 ) * ( . 5 * RHO -~L5 / S ) 

YDDD  = ( -  0 . 0005 ) * ( . 5 -RHO *L2 * S * * 2 ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

C      YVDOT=.(  -0  .  0039  )*(  .  5-RHO-L3  ) 
C   SPEED=23  KNOTS,  ENCOUNTER  FREQUENCY  =0.75 

YVDOT=-2304300.0 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV  = ( -  0 . 0  0  2 1 3 ) * ( . 5 -RHO  * L  3 * S ) 

NR= ( -  0 . 00 10  5 ) * ( . 5 -RHO *L4 * S ) 

ND= ( -  0 . 0007 ) - ( . 5 -RHO-L3 - S - -2 ) 

NVVR  = ( -  0 . 0 15 ) * ( . 5 -RHO -L4 / S ) 

NVRR= (-0.008) * ( . 5 * RHO * L5 / S ) 
' NVVV= ( 0 . 0 1 ) * ( . 5 -RHO -L3 / S ) 

NRRR= ( -  0 . 0  0  6 ) - ( . 5  * RHO * L 6 / S ) 

NDDD  = ( 0 . 000 1 ) - ( . 5 -RHO -L3 - S - -2 ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
C 

C      NRDOT  = (-0.00027)* (.5 * RHO *L5 ) 
C   SPEED=23  KNOTS,  ENCOUNTER  FREQUENCY  =0.75 

NRDOT=-1.4518E+ll 
C   SETS  SEA  STATE  TO  ZERO 

IF  (ISEA.EQ. 1)  GO  TO  30 

FX=0. 

FY=0. 

MZ  =  0. 
GO  TO  35 
C   UNIT  12  HAS  THE  SEA  STATE  DATA  NAMED  CH 
C   IT  MUST  BE  SYNCHRONIZED  BY  TIME 
30    READ  (12)  CH 

FX=CH(3) 

FY=CH(4) 

MZ=CH(8) 
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35    CONTINUE 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU*UC**2 
C   EQUATIONS  OF  MOTION 

C      UDOT=(  (XVR  +  MASS)*V*R  +  XUU*U**2  +  XW*V**2 
C     1  +  XDD*D*D  +   FX  +  XP  ) / (MASS-XUDOT ) 

VDOT=(YV*V  +  (YR-MASS*U)*R  +  YD*D  +  YVVR*V**2*R 

1  +  YVRR*V*R**2  +  YVW*V**3 

2  +  YRRR*R**3  +  YDDD*D**3  +  FY  )  /  (MASS  -  YVDOT ) 
RDOT=(NV*V  +  NR-'R  +  ND-D  +  NVVR*V**2*R  +  NVRR*V*R**2 

I  +  NVW*V**3  +  NRRR*R**3  +  NDDD*D**3  +  MZ )  /  ( IZ-NRDOT ) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 

GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
5  0   YAWDEG  =  YAW- 5  7.296 

RDEG=R*5  7.29  6 

RDDEG=RDOT*57.296 

DDEG=D*57.296 

YAWC=YAWC*57.296 

ICOUNT=l 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1.0 
C   INTEGRATION 

U=U+UDOT*DELT 

V=V+VDOT*DELT 

R=R+RDOT-DELT 

YAW=YAW+R*DELT 

X2=X2+DX2*DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C      XDOT=U"DCOS(YAW)-V-DSIN(YAW) 
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C      YDOT  =  U-DS IN ( YAW ) + V -DCOS ( YAW ) 
C      X=X+XDOT*DELT 
C      Y=Y+YDOT-DELT 

TIME=TIME+DELT 

ICOUNT=ICOUNT+l 

ISE=ISE  +  LAMDA*YAWE**2 

ISR=ISR  +  D**2 

GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE(6,500)  TDIFF,K1,T1,T2 
500   FORMAT ('  *, IX, 'TOTAL  =',F15.7,?   Kl  =',F15.7, 
1  *  Tl  = ' ,F15.7,2X, 'T2=' ,F15.7) 

REWIND  12 

RETURN 

END 
The  function  minimization  subroutine  BOXPLX  follows 
Then  the  following  three  cards  must  be  placed. 
//GO. SYSIN  DD  * 

/* 

/  GO.FT12F001  DD  DISP=SHR, DSN=MSS . S2 160 . A241 
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APPENDIX  E 
SYSTEM'S  RESPONSE  FOR  IRREGULAR  SEAS 

//IRRERESP  JOB  (XXXX , XXXX) ,' RESEARCH ', CLASS^B 

//-MAIN  ORG=NPGVMl.XXXXP 

//  EXEC  FRTXCLGP,IMSL=DP,REGION=1024K 

//FORT. SYSIN  DD  * 

C   IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 

C   OBTAINED 

DIMENSION  XX(3) 
C   OPTIMAL  GAINS  FOR  CONTROLLER 
XX(1)=. 9192610 
XX(2)=18. 5788116 
XX(3)=9. 77668 
C   THE  SUBROUTINE  PLANT  SIMULATES  THE  SL-7  CONTAINERSHIP 

CALL  PLANT (XX) 

WRITE (6, 25) 
25     FORMAT (IX, 'OPTIMAL  GAINS' ,/ ) 

DO  30  1=1,3 
30     WRITE (6 ,40) I, XX (I) 
40     F0RMAT(1X, 'XX( ' ,12, ' )=' ,F14.7) 

STOP 

END 
C 
C   SUBROUTINE  PLANT (XX)   SIMULATES  THE  SHIP 

SUBROUTINE  PLANT (XX) 

COMMON  TDIFF 

REAL-8  L,L2,L3,L4,L5 ,L6 

REAL-8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 

REAL- 8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 

REAL- 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 

REAL -8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL- 8  RHO , IZ , FX , FY , MZ , XP , MASS , DELT 
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REAL- 8  DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D 

REAL- 8  K1,T1,T2,D,X2,DX2,S,CH(11) ,DX3 ,X3,X4 

DIMENSION  XX(3) 
C 

C   CLOSE  LOOP  ANALYSIS  WITH  FILTER 
C 

INITIAL  CONDITIONS  FOR  INTEGRATION 
C   SIMULATION  END  TIME  IN  SECONDS 

ETIME=600. 

TIME=0.0 

ICOUNT=l 
C   INITIALIZE  THE  COST  FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=4.2 
C   GAIN  COEFFICIENTS 

K1=XX(1) 

T1=XX(2) 

T2=XX(3) 
C   X,XDOT,Y,YDOT  ARE  FIX  COORDINATES  ON  EARTH 

X=0.0 

Y=0.0 

XDOT=0.0 

YDOT=0.0 
C   U,UDOT,V,VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V=0.0 

UDOT=0.0 

VDOT=0.0 

YAW^O.O 

R=0.0 

RDOT=0.0 

YAW^O.O 
C   ORDERED  SPEED  IN   FEET/ SEC 
C    38.81  FT/SEC=23  KNOTS 
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UC-38.81 
C   AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 
C   D  =  RUDDER  ANGLE 

D=0.0 

L=880.5 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 
C   SEA  DISTURBANCE 

C   FORCES  IN  XJ  DIRECTION.  COMPUTED  IN  FORCES 
C   MOMENTS  IN  Z 

FX=0. 

FY=0. 

MZ  =  0. 
C   ISEA  IS  A  SWITCH;  ISEA=0  (CALM  WATER)  ISEA=1  (SEA  STATE) 

ISEA=1 
C   HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 

RHO=1.9876 

MASS=( .0044)*( .5*RHO*L3) 

IZ=(0.00028)*(.5 * RHO - L5 ) 

YAWE=0.0 

X2=0.0 

DX2=0.0 

X3=0.0 

DX3=0.0 

X4=0.0 
200  CONTINUE 

S=DSQRT(U**2  +  V**2) 
C   INPUT  YAW  COMMAND 

YAWC=0.0 

IF  (TIME. GE. 0.0)  YAWC=0.0 
C   ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 
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C   (  COMPENSATOR  FILTER  ) 
YAWE=YAW  -  YAWC 
DX2=(YAWE-X2)/T2 
X4=K1*(T1"DX2+X2) 

DX3=(X4-X3)/T4 

D=(T3-DX3+X3) 
C   AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
C 

C   XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 
C      XUDOT  = ( - . 000 1 ) * ( . 5 * RHO *L3  ) 

XUU= ( -  0 . 0  0  0  3 ) * ( . 5 * RHO *L2 ) 

XVR=  (0.0039  ).*  (  .  5  *  RHO  -  L3  ) 

XVV= ( - . 00 12 ) * ( . 5 -'RHO -'L2 ) 

XDD= ( -  0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 
C   LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

YV= ( -0 . 00758 )*( . 5*RHO*L2*S ) 

YR= ( 0 . 0  0  2  3 ) * ( . 5 * RHO - L3  * S ) 

YD= ( 0 . 00145 ) * ( . 5 "RH0-L2-S * *2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 008 ) * ( . 5 -RHO -L4 / S ) 

YWV=  (  -  0  .  0 3  )  -  (  .  5  -RHO -L2 /  S  ) 

YRRR= ( 0 . 0  0  3 ) * ( . 5 -RHO * L5 / S ) 

YDDD= ( -  0 . 0005 ) * ( . 5 -RHO-L2 * S * -2 ) 
C   YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C   DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 
C 
C      YVDOT  = ( -  0 . 0  0  3  9 ) * ( . 5 * RHO * L3 ) 

YVDOT=-3654800.00 
C   MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= (  -  0 . 0 0 2 1 3 ) - ( . 5  - RHO - L3 - S  ) 

NR= ( -  0 . 0  0 1 0  5 ) * ( . 5  - RHO - L4 - S ) 

ND  =  (  -  0  .  0  0  0  7 ) - ( . 5  - RHO * L 3  * S  * * 2  ) 

NVVR= ( -  0 . 0 15 ) - ( . 5 -RHO -L4/ S ) 

NVRR= ( -  0 . 0 0 8 ) - ( . 5  - RHO - L5 / S ) 
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NVVV  =  (  0  .  0 1 )  *  (  .  5  * RHO  * L  3  /  S  ) 
NRRR= ( -  0 . 0  0  6 ) * ( . 5 -RHO *L 6 / S ) 
NDDD  = ( 0 . 000 1 ) * ( . 5*RHO*L3*S**2 ) 
C   NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C   FOR  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 
C 
C      NRDOT  = ( -  0 . 0  0  0 2  7 ) * ( . 5 * RHO * L5 ) 

NRDOT=-2. 1815E+11 
C   SETS  SEA  STATE  TO  ZERO 

IF  (ISEA.EQ. 1)  GO  TO  30 
FX=0. 
FY=0. 
MZ  =  0. 
GO  TO  35 
C   UNIT  12  HAS  THE  SEA  STATE  DATA  NAMED  CH 
C   IT  MUST  BE  SYNCHRONIZED  BY  TIME 
30     READ  (12)  CH 
FX=  CH(3) 
FY=  CH(4) 
MZ  =  CH(8) 
35     CONTINUE 
C   U  ACTUAL  SPEED 
C   UC  COMMANDED  SPEED 
C   XP  =  PROPELLER  THRUST 

XP=-XUU*UC**2 
C   EQUATIONS  OF  MOTION 

C      UDOT=(  (XVR  +  MASS)*V*R  +  XUU*U**2  +  XW*V**2 
C     1  +  XDD-D-D  +  FX   +  XP  ) / (MASS-XUDOT ) 

VDOT=(YV*V  +  (YR-MASS*S)*R  +  YD-D  +  YVVR*V**2*R 
1  +  YVRR*V*R**2  +  YWV*V**3 
1  +  YRRR*R**3  +  YDDD-D--3  +  FY  ) / (MASS  - YVDOT ) 

RDOT=(NV*V  +  NR"R  +  ND--D  +  NWR*V**2*R  +  NVRR*V*R**2 
1  +  NVVV*V**3  +  NRRR*R**3  +  NDDD-D--3  +  MZ )/( IZ- NRDOT ) 
C   WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.2  )  GO  TO  50 
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GO  TO  300 
C   CONVERT  RADIANS  TO  DEGREES 
5  0   YAWDEG  =  YAW * 5  7.296 

RDEG=R*57.296 

RDDEG=RDOT*57.296 

DDEG=D*57.296 

YAWC=YAWC-5  7.29  6 

WRITE  (6,100)  TIME, YAWDEG 
100  FORMAT ( IX, F 12. 8 ,1X,F12.8) 

ICOUNT=l 
C   TEST  IF  WANT  TO  STOP 

300   IF  (TIME.GE.ETIME)  GO  TO  400 
C   INTEGRATION  STEP  SIZE  DELT 

DELT=1. 
C   INTEGRATION 

U=U+UDOT*DELT 

V=V+VDOT*DELT 

R=R+RDOT*DELT 

YAW=YAW+R*DELT 

X2=X2+DX2*DELT 

X3=X3+DX3*DELT 
C   CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 

XDOT=U»DCOS(YAW)-V-DSIN(YAW) 

YDOT=U*DSIN(YAW)+V*DCOS(YAW) 

X=X+XDOT*DELT 

Y=Y+YDOT*DELT 

TIME=TIME+DELT 

ICOUNT=ICOUNT+l 

ISE  =  ISE  +  LAMDA*YAWE**2 

ISR=ISR  +  D**2 

GO  TO  200 
C   J=TDIFF=  COST  FUNCTION 
400   TDIFF=ISE+ISR 

WRITE(6,500)  ISE,ISR,TDIFF 
500   FORMAT( '1'  ,5X,  ' ISE= '  ,F15.7,  '   ISR=',F15.7, 
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1  '  TOTALS  '  ,F15.7) 
STOP 
END 
//GO. SYS IN  DD  * 

/* 

//GO.FT12F001  DD  DISP=SHR ,DSN=MSS . S2160 . A211 
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APPENDIX  F 
MODIFIED  MINIMIZATION  SUBROUTINE 

SUBROUTINE  BOXPLX  (NV , NAV , NPR , NTZ , RZ , XS , IP , BU , BL , YMN , IER) 
C 

DIMENSION  V(50, 50) ,  FUN(50),  SUM(25),  CEN(25),  XS(NV), 
1BU(NV) ,BL(NV) 
C 

KV  =    5 

EP  =  l.E-4 

NTA  =  2000 

IF  (NTZ.GT.O)  NTA  =  NTZ 

R  =  RZ 

IF  (R.LE.O. .OR.R.GE. 1. )  R= 1 . / 3 . 

NVT  =  NV+NAV 
C 
C     TOTAL  VARS,  EXPLICIT  PLUS  IMPLICIT 

NT  =  0 
C     CURRENT  TRIAL  NO. 

NPT  =  0 
C     CURRENT  NO.  OF  PERMISSIBLE  TRIALS 

NTFS  =    0 
C     CURRENT  NO.  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGED 
C 
C     CHECK  FEASIBILITY  OF  START  POINT 


C 


DO  4  I=1,NV 

VT  =  XS(I) 

IF  (BL(I) .LE.VT)  GO  TO  1 

II  =  -I 

VT  =  BL(I) 

GO  TO  2 

1  IF  (BU(I) .GE.VT)  GO  TO  3 
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c 
c 


II  =  I 

VT  =  BU(I) 

2  IF  (NPR.GT.O)  WRITE  (6,49)  II 

3  V(I,1)  =  VT 
CEN(I)  =  VT 

IF  (IP.EQ. 1)  GO  TO  4 

BL(I)  =  BL(I)+AMAX1(EP,EP*ABS(BL(I))) 

BU(I)  =  BU(I)-AMAX1(EP,EP*ABS(BU(I))) 

4  SUM(I)  =  VT 


NCE  =  1 
C     NUMBER  OF  CONSTRAINT  EVALUATIONS 

I  =  1 

IF  (KE(V(1,1) ) .EQ.O)  GO  TO  5 

IF  (NPR.LE.O)  GO  TO  12 

WRITE  (6,50) 

GO  TO  12 
5  NFE  =  1 
C 
C    NUMBER  OF  VERTICES  (K)  =  2  TIMES  NO.  OF  VARIABLES 

K  =  (2*NV)/3 
C 
C    NUMBER  OF  DISPLACEMENTS  ALLOWED. 

NLIM  =  5-NV+10 
C 

C    NUMBER  OF  CONSECUTIVE  TRIALS  WITH  UNCHANGED  FE  TO 
C    TO  TERMINATE 

NCT  =  NLIM+NV 

ALPHA  =1.3 

FK  =  K 

FKM  =  FK-1. 

BETA  =  ALPHA+1. 
C 
C    INSURE  SEED  OF  RANDOM  NUMBER  GENERATOR  IS  ODD. 

103 


IQR  =  R-1.E7 

IF  (M0D(IQR,2) .EQ.O)  IQR=IQR+101 
C 
C     SET  UP  INITIAL  VERTICES 

FUN(l)  =  FE(V(1,1)) 

YMN  =  FUN(l) 
6  FI  =  1. 

FUNOLD  =  FUN(l) 


DO  15  1  =  2, K 

FI  =  FI+1. 

LIMT  =  0 

7  LIMT  =  LIMT+1 

c 

c 

END  CALCULATIO 

IF  (LIMT.GE.NLIM)  GO  TO  11 
C 

DO  8  J=1,NV 
C 
C    RANDOM  NUMBER  GENERATOR   (RANDU) 

IQR  =  IQR- 6553 9 

IF  (IQR.LT.O)  IQR  =  IQR+2147483647+ 1 

RQX  =  IQR 

RQX  =  RQX*.4656613E-9 

V(J,I)  =  BL(J)+RQX*(BU(J)-BL(J) ) 

IF  (IP.EQ.l)  V(J,I)=AINT(V(J,I)+.5) 
8  CONTINUE 
C 

DO  10  L=1,NLIM 

NCE  =  NCE+1 

IF  (KE(V(1,I)) .EQ.O)  GO  TO  13 
C 

DO  9  J=1,NV 

VT  =  .5*(V(J,I)+CEN(J)) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 
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V(J,I)  =  VT 
9  CONTINUE 
C 

10  CONTINUE 
C 

11  IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,51)  I 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,I,FUN,CEN,I) 

12  IER  =  -1 
GO  TO  48 

C 

13  DO  14  J=1,NV 

SUM(J)  =  SUM(J)+V(J,I) 

14  CEN(J)  =  SUM(J)/FI 
C 

C    TRY  TO  ASSURE  FEASIBLE  CENTROID  FOR  STARTING. 
NCE  =  NCE+1 

IF  (KE(CEN) .EQ.O)  GO  TO  60 
SUM(J)  =  SUM(J)  -V(J,I) 
GO  TO  7 
60  NFE  =  NFE+1 

FUN(I)  =  FE(V(1,I)) 

15  CONTINUE 
C 

C    END  OF  LOOP  SETTING  OF  INITIAL  COMPLEX. 

IF  (NPR.LE.O)  GO  TO  17 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 
C 
C    FIND  THE  WORST  VERTEX,  THE  ' J ' TH . 

J  =  1 
C 

DO  16  1=2, K 

IF  (FUN(J) .GE.FUN(I) )  GO  TO  16 

J  =  I 

16  CONTINUE 
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c 

C    BASIC  LOOP.   ELIMINATE  EACH  WORST  VERTEX  IN  TURN. 
C    IT  MUST  BECOME  NO  LONGER  WORST,  NOT  MERELY  IMPROVED 
C    FIND  NEXT-TO-WORST  VERTEX,   THE  ' JN ' TH  ONE. 

17  JN  =  1 

IF  (J.EQ. 1)  JN  =  2 
C 

DO  18  1=1, K 

IF  (I.EQ.J)  GO  TO  18 

IF  (FUN(JN) .GE.FUN(I) )  GO  TO  18 

JN  =  I 

18  CONTINUE 
C 

C    LIMT  =  NUMBER  OF  MOVES  DURING  THIS  TRIAL  TOWARD  THE 
C    CENTROID  DUE  TO  FUNCTION  VALUE. 

LIMT  =  1 
C 

C    COMPUTE  CENTROID  AND  OVER  REFLECT  WORST  VERTEX. 
C 

DO  19  1=1, NV 

VT  =  V(I,J) 

SUM(I)  =  SUM(I)-VT 

CEN(I)  =  SUM(I)/FKM 

VT  =  BETA*CEN(I)-ALPHA*VT 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 
C 
C    INSURE  THE  EXPLICIT  CONSTRAINTS  ARE  OBSERVED. 

19  V(I,J)  =  AMAX1(AMIN1(VT,BU(I) ) ,BL(I) ) 
C 

NT  =  NT+1 
C 
C    CHECK  FOR  IMPLICIT  CONSTRAINT  VIOLATION. 


C 


20  DO  25  N=1,NLIM 
NCE  =  NCE+1 
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IF  (KE(V(1,J)) .EQ.O)  GO  TO  26 
C 

C    EVERY  'KV'TH  TIME,  OVER-REFLECT  THE  OFFENDING  VERTEX 
C    THROUGH  THE  BEST  VERTEX. 

IF  (MOD(N,KV) .NE.O)  GO  TO  22 

CALL  FBV  (K,FUN,M) 
C 

DO  21  1=1, NV 

VT  =  BETA"V(I,M)-ALPHA-V(I,J) 

IF  (IP.EQ.l)  VT  =    AINT(VT+.5) 

21  V(I,J)  =    AMAX1(AMIN1(VT,BU(I) ) ,BL(I)) 
C 

GO  TO  24 
C 

C    CONSTRAINT  VIOLATION:   MOVE  NEW  POINT  TOWARD  CENTROID. 
C 

22  DO  23  1=1, NV 

VT  =  . 5*(CEN(I)+V(I, J) ) 

IF  (IP.EQ.l)  VT  =  AIN.T(VT+.5) 

V(I,J)  =  VT 

23  CONTINUE 
C 

24  NT  =  NT+1 

25  CONTINUE 
C 

IER  =  1 
C 

C    CANNOT  GET  FEASIBLE  VERTEX  BY  MOVING  TOWARD  CENTROID, 
C    OR  BY  OVER-REFLECTING  THRU  THE  BEST  VERTEX. 

IF  (NPR.LE.O)  GO  TO  42 

WRITE  (6,52)  NT, J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 

GO  TO  42 
C 
C    FEASIBLE  VERTEX  FOUND,  EVALUATE  THE  OBJECTIVE  FUNCTION 
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26  NFE  =  NFE+1 
FUNTRY  =  FE(V(1,J)) 

C 

C    TEST  TO  -  SEE  IF  FUNCTION  VALUE  HAS  NOT  CHANGED. 

AFO  =  ABS( FUNTRY- FUNOLD) 

AMX  =  AMAX  1(ABS(EP-- FUNOLD )  ,EP) 
C 

C    ACTIVATE  THE  FOLLOWING  TWO  STATEMENTS 
C    FOR  DIAGNOSTIC  PURPOSES  ONLY. 

C      WRITE  (6,99)  J , AFO , AMX , FUNTRY , FUNOLD , FUN (J ), FUN (JN 
C     1NTFS,N 
C   99  FORMAT  ( IX , 13 , 6E15 . 7 , 215 ) 

IF  (AFO. GT. AMX)  GO  TO  27 

NTFS  =  NTFS+1 

IF  (NTFS.LT.NCT)  GO  TO  28 

IER  =  0 

IF  (NPR.LE.O)  GO  TO  42 

WRITE  (6,53)  K 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 

GO  TO  42 

27  NTFS  =  0 
C 

C    IS  THE  NEW  VERTEX  NO  LONGER  WORST? 

28  IF  ( FUNTRY. LT. FUN (JN))  GO  TO  34 
C 

C  TRIAL  VERTEX  IS  STILL  WORST;  ADJUST  TOWARD  CENTROID. 
C  EVERY  'KV'TH  TIME,  OVER- REFLECT  THE  OFFENDING  VERTEX 
C    THROUGH  THE  BEST  VERTEX. 

LIMT  =  LIMT+1 

IF  (MOD (LIMT, KV ) . NE . 0 )  GO  TO  30 

CALL  FBV  (K,FUN,M) 
C 

DO  29  1=1, NV 

VT  =  BETA*V(I,M)-ALPHA*V(I,J) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 
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29  V(I,J)  =  AMAX1(AMIN1(VT,BU(I) ) ,BL(I)) 
C 

GO  TO  32 
C 

30  DO  31  1=1, NV 

VT  =  .5*(CEN(I)+V(I,J)) 

IF  (IP.EQ.l)  VT  =  AINT(VT+.5) 

V(I,J)  =  VT 

31  CONTINUE 
C 

32  IF  (LIMT.LT.NLIM)  GO  TO  33 
C 

C    CANNOT  MAKE  THE  ' J ' TH  VERTEX  NO  LONGER  WORST  BY 

C    DISPLACING  TOWARD  THE  CENTROID  OR  BY  OVER-REFLECTING 

C    THRU  THE  BEST  VERTEX. 

IER  =  2 

IF  (NPR  . LE.  0)   GO  TO  42 

WRITE  (6,52)   NT,  J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 

GO  TO  42 

33  NT  =  NT+1   • 
GO  TO  20 

C 

C    SUCCESS:   WE  HAVE  A  REPLACEMENT  FOR  VERTEX  J. 

34  FUN (J)  =  FUNTRY 
FUNOLD  =  FUNTRY 
NPT  =  NPT+1 

C 

C    STOP  AT  THE  100' TH  PERMISSIBLE  TRIAL 

IF  (MOD(NPT, 100) .EQ.O)  GO  TO  48 
C 

DO  3  6  1=1, NV 

SUM(I)  =  0. 
C 

DO  35  N=1,K 
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3  5  SUM(I)  =  SUM(I)+V(I,N) 
C  • 

CEN(I)  =  SUM(I)/FK 

36  CONTINUE 
C 

LC  =  0 
GO  TO  39 
C 

37  DO  38  1=1, NV 

38  SUM(I)  =  SUM(I)+V(I,J) 
C 

LC  =  J 
C 

39  IF  (NPR.LE.O)  GO  TO  40 

IF  (MOD(NPT,NPR) . NE . 0 )  GO  TO  40 
C 

CALL  BOUT  (NT,N?T,NFE,NCE,NV,NVT,V,K,FUN,CEN,LC) 
C 

C    HAS  THE  MAX.  NUMBER  OF  TRIALS  BEEN  REACHED  WITHOUT 
C    CONVERGENCE  ?  IF  NOT,  GO  TO  NEW  TRIAL. 

40  IF  (NT.GE.NTA)  GO  TO  41 
C 

C        NEXT-TO-WORST  VERTEX  NOW  BECOMES  WORST. 
J  =  JN 
GO  TO  17 

41  IER  =  3 

IF  (NPR.GT.O)  WRITE  (6,54) 
C 

C    COLLECTOR  POINT  FOR  ALL  ENDINGS. 

C     1)   CANNOT  DEVELOP  FEASIBLE  VERTEX.  IER  =  1 

C  2)  CANNOT  DEVELOP  A  NO-LONGER- WORST  VERTEX.  IER  =  2 
C  3)  FUNCTION  VALUE  UNCHANGED  FOR  K  TRIALS.  IER  =  0 
C     4)   LIMIT  ON  TRIALS  REACHED.  IER  =  3 

C     5)   CANNOT  FIND  FEASIBLE  VERTEX  AT  START.      IER  =-1 

42  CONTINUE 
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c 

C    FIND  BEST  VERTEX. 

CALL  FBV  (K,FUN,M) 

IF  (IER.GE.3)  GO  TO  44 
C    RESTART  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER  THAN 
C    THE  PREVIOUS  OR  IF  THIS  IS  THE  FIRST  TRY. 

IF  (NPR.LE.O)  GO  TO  43 

WRITE  (6,55)  (M,YMN,FUN(M)) 

43  IF  (FUN(M) .GE.YMN)  GO  TO  47 

IF  (ABS(FUN(M)-YMN)  .  LE  .  AMAX1  (EP  ,  EP-'YMN)  )  GO  TO  47 
C    GIVE  IT  ANOTHER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 

44  YMN  =    FUN(M) 
FUN(l)  =  FUN(M) 

C 

DO  45  1=1, NV 
CEN(I)  =  V(I,M) 
SUM(I)  =    V(I,M) 

45  V(I,1)  =  V(I,M) 
C 

DO  46  1=1, NVT 

46  XS(I)  =  V(I,M) 
C 

IF  (IER.LT.3)  GO  TO  6 

47  IF  (NPR.LE.O)  GO  TO  48 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,V(1,M) , -1) 
WRITE  (6,56)  FUN(M) 

48  RETURN 
C 

49  FORMAT  (50H0INDEX  AND  DIRECTION  OF  OUTLYING 
1VARIABLE  AT  STARTI5) 

50  FORMAT  (50H0IMPLICIT  CONSTRAINT  VIOLATED  AT  START. 
1DEAD  END . ) 

51  FORMAT  ('OCANNOT  FIND  FEASIBLE ', 14 ,' TH  VERTEX  OR 
1CENTROID  AT  START.') 

52  FORMAT  (10H0AT  TRIAL  I4,54H  CANNOT  FIND  FEASIBLE 
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1VERTEX  WHICH  IS  NO 

2LONGER  WORST, 14, 15X, 'RESTART  FROM  BEST  VERTEX.') 
5  3  FORMAT  (40H0FUNCTION  HAS  BEEN  ALMOST  UNCHANGED 
IFOR  15 ,7H  TRIALS) 

54  FORMAT  (27H0LIMIT  ON  TRIALS  EXCEEDED.  ) 

55  FORMAT  ( ' OBEST  VERTEX  IS  NO. ',13,'   OLD  MIN  WAS  ', 
1E15.7,  '   NEW  MIN  IS  ',E15.7) 

56  FORMAT  ( ' OMIN  OBJECTIVE  FUNCTION  IS  ',E15.7) 
END 

SUBROUTINE  FBV  (K,FUN,M) 
DIMENSION  FUN (50) 
M  =  1 

DO  1  1  =  2, K 

IF  (FUN(M) .LE.FUN(I) )  GO  TO  1 
M  =  I 
1  CONTINUE 

RETURN 

END 

SUBROUTINE  BOUT  (NT , NPT ,NFE , NCE , NV , NVT , V ,K , FN , C , IK) 

DIMENSION  V(50, 50) ,  FN(50),  C(25) 

WRITE  (6,4)  NT,NPT,NFE,NCE 

DO  1  1=1, K 

WRITE  (6,5)  FN(I) , (V(J,I) ,J=1,NV) 
IF  (NVT.LE.NV)  GO  TO  1 
NVP  =  NV+1 

WRITE  (6,6)  (V(J,I) ,J=NVP,NVT) 
1  CONTINUE 

IF  (IK.NE.O)  GO  TO  2 

WRITE  (6,7)  (C(I),I=1,NV) 

RETURN 
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2  IF  (IK.GE.O)  GO  TO  3 
WRITE  (6,8)  (C(I) ,I=1,NV) 
RETURN 

3  WRITE  (6,9)  IK, (C(I) ,I=1,NV) 
RETURN 

4  FORMAT  ('ONO.  TOTAL  TRIALS  =  ',I5,4X,'NO.  FEASIBLE 
1TRIALS  =  •  '  ,' 

2I5,4X,'N0.  FUNCTION  EVALUATIONS  =  ',I5,4X,'NO. 

3CONSTRAINT  EVALUATIONS 

4=  ',I5/'0      FUNCTION  VALUE' ,6X, 'INDEPENDENT 

5VARIABLES / DEPENDENT 

60R  IMPLICIT  CONSTRAINTS') 

5  FORMAT  (1H  , E18 . 7 , 2X , 7E14 . 7/ ( 21X , 7E14 . 7 ) ) 

6  FORMAT  (21X,7E14.7) 

7  FORMAT  (10H0CENTROID  11X , 7E14 . 7 / (21X , 7E14 . 7 ) ) 

8  FORMAT  ('0   BEST  VERTEX ' , 7X , 7E14 . 7/ (21X , 7E14 . 7 ) ) 

9  FORMAT  ('OCENTROID  LESS  VX ' , 12 , 2X , 7E14 . 7/ ( 2IX , 7E14 . 7 ) ) 
END 

FUNCTION  FE(X) 

DIMENSION  X(3)- 

COMMON  TDIFF 

CALL  PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION  KE(X) 

DIMENSION  X(3) 

KE  =  0 

RETURN 

END 
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